(git:ccc2433)
qs_moments.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 !> \brief Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
10 !> \par History
11 !> added angular moments (JGH 11.2012)
12 !> \author JGH (20.07.2006)
13 ! **************************************************************************************************
14 MODULE qs_moments
15  USE ai_angmom, ONLY: angmom
16  USE ai_moments, ONLY: contract_cossin,&
17  cossin,&
18  diff_momop,&
19  diff_momop2,&
20  diff_momop_velocity,&
21  moment
22  USE atomic_kind_types, ONLY: atomic_kind_type,&
23  get_atomic_kind
24  USE basis_set_types, ONLY: gto_basis_set_p_type,&
25  gto_basis_set_type
26  USE bibliography, ONLY: mattiat2019,&
27  cite_reference
28  USE block_p_types, ONLY: block_p_type
29  USE cell_types, ONLY: cell_type,&
30  pbc
31  USE commutator_rpnl, ONLY: build_com_mom_nl
32  USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
33  USE cp_cfm_types, ONLY: cp_cfm_create,&
34  cp_cfm_get_info,&
35  cp_cfm_release,&
36  cp_cfm_type
37  USE cp_control_types, ONLY: dft_control_type
38  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
39  USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
40  dbcsr_allocate_matrix_set,&
41  dbcsr_deallocate_matrix_set
42  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43  cp_fm_struct_double,&
44  cp_fm_struct_release,&
45  cp_fm_struct_type
46  USE cp_fm_types, ONLY: cp_fm_create,&
47  cp_fm_get_info,&
48  cp_fm_release,&
49  cp_fm_set_all,&
50  cp_fm_type
51  USE cp_result_methods, ONLY: cp_results_erase,&
52  put_results
53  USE cp_result_types, ONLY: cp_result_type
54  USE dbcsr_api, ONLY: &
55  dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_dot, &
56  dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_trace, &
57  dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
58  USE distribution_1d_types, ONLY: distribution_1d_type
59  USE kinds, ONLY: default_string_length,&
60  dp
61  USE kpoint_types, ONLY: get_kpoint_info,&
62  kpoint_type
63  USE mathconstants, ONLY: pi,&
64  twopi
65  USE message_passing, ONLY: mp_para_env_type
66  USE moments_utils, ONLY: get_reference_point
67  USE orbital_pointers, ONLY: current_maxl,&
68  indco,&
69  ncoset
70  USE parallel_gemm_api, ONLY: parallel_gemm
71  USE particle_methods, ONLY: get_particle_set
72  USE particle_types, ONLY: particle_type
73  USE physcon, ONLY: bohr,&
74  debye
75  USE qs_environment_types, ONLY: get_qs_env,&
76  qs_environment_type
77  USE qs_kind_types, ONLY: get_qs_kind,&
78  get_qs_kind_set,&
79  qs_kind_type
80  USE qs_ks_types, ONLY: get_ks_env,&
81  qs_ks_env_type
82  USE qs_mo_types, ONLY: get_mo_set,&
83  mo_set_type
84  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
85  neighbor_list_iterate,&
86  neighbor_list_iterator_create,&
87  neighbor_list_iterator_p_type,&
88  neighbor_list_iterator_release,&
89  neighbor_list_set_p_type
90  USE qs_operators_ao, ONLY: build_lin_mom_matrix
91  USE qs_rho_types, ONLY: qs_rho_get,&
92  qs_rho_type
93  USE rt_propagation_types, ONLY: get_rtp,&
94  rt_prop_type
95 #include "./base/base_uses.f90"
96 
97  IMPLICIT NONE
98 
99  PRIVATE
100 
101  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
102 
103  ! Public subroutines
104  PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
105  PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
106  PUBLIC :: qs_moment_berry_phase, qs_moment_locop
107  PUBLIC :: dipole_deriv_ao
108  PUBLIC :: build_local_moments_der_matrix
109  PUBLIC :: build_dsdv_moments
110  PUBLIC :: dipole_velocity_deriv
111 
112 CONTAINS
113 
114 ! *****************************************************************************
115 !> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
116 !> to the basis function on the right
117 !> difdip(beta, alpha) = < mu | r_beta | ∂_alpha nu > * (mu - nu)
118 !> \param qs_env ...
119 !> \param difdip ...
120 !> \param order The order of the derivative (1 for dipole moment)
121 !> \param lambda The atom on which we take the derivative
122 !> \param rc ...
123 !> \author Edward Ditler
124 ! **************************************************************************************************
125  SUBROUTINE dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
126  TYPE(qs_environment_type), POINTER :: qs_env
127  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: difdip
128  INTEGER, INTENT(IN) :: order, lambda
129  REAL(kind=dp), DIMENSION(3) :: rc
130 
131  CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_velocity_deriv'
132 
133  INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
134  last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
135  sgfb
136  LOGICAL :: found
137  REAL(dp) :: dab
138  REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
139  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
140  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
141  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab, difmab2
142  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
143  TYPE(cell_type), POINTER :: cell
144  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
145  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
146  TYPE(neighbor_list_iterator_p_type), &
147  DIMENSION(:), POINTER :: nl_iterator
148  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149  POINTER :: sab_all
150  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
152  TYPE(qs_kind_type), POINTER :: qs_kind
153 
154  CALL timeset(routinen, handle)
155 
156  NULLIFY (cell, particle_set, qs_kind_set, sab_all)
157  CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
158  qs_kind_set=qs_kind_set, sab_all=sab_all)
159  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
160  maxco=ldab, maxsgf=maxsgf)
161 
162  nkind = SIZE(qs_kind_set)
163  natom = SIZE(particle_set)
164 
165  m_dim = ncoset(order) - 1
166 
167  ALLOCATE (basis_set_list(nkind))
168 
169  ALLOCATE (mab(ldab, ldab, m_dim))
170  ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
171  ALLOCATE (work(ldab, maxsgf))
172  ALLOCATE (mint(3, 3))
173  ALLOCATE (mint2(3, 3))
174 
175  mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
176  difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
177  work(1:ldab, 1:maxsgf) = 0.0_dp
178 
179  DO i = 1, 3
180  DO j = 1, 3
181  NULLIFY (mint(i, j)%block)
182  NULLIFY (mint2(i, j)%block)
183  END DO
184  END DO
185 
186  ! Set the basis_set_list(nkind) to point to the corresponding basis sets
187  DO ikind = 1, nkind
188  qs_kind => qs_kind_set(ikind)
189  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
190  IF (ASSOCIATED(basis_set_a)) THEN
191  basis_set_list(ikind)%gto_basis_set => basis_set_a
192  ELSE
193  NULLIFY (basis_set_list(ikind)%gto_basis_set)
194  END IF
195  END DO
196 
197  CALL neighbor_list_iterator_create(nl_iterator, sab_all)
198  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
199  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
200  iatom=iatom, jatom=jatom, r=rab)
201 
202  basis_set_a => basis_set_list(ikind)%gto_basis_set
203  basis_set_b => basis_set_list(jkind)%gto_basis_set
204  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
205  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
206 
207  associate( &
208  ! basis ikind
209  first_sgfa => basis_set_a%first_sgf, &
210  la_max => basis_set_a%lmax, &
211  la_min => basis_set_a%lmin, &
212  npgfa => basis_set_a%npgf, &
213  nsgfa => basis_set_a%nsgf_set, &
214  rpgfa => basis_set_a%pgf_radius, &
215  set_radius_a => basis_set_a%set_radius, &
216  sphi_a => basis_set_a%sphi, &
217  zeta => basis_set_a%zet, &
218  ! basis jkind, &
219  first_sgfb => basis_set_b%first_sgf, &
220  lb_max => basis_set_b%lmax, &
221  lb_min => basis_set_b%lmin, &
222  npgfb => basis_set_b%npgf, &
223  nsgfb => basis_set_b%nsgf_set, &
224  rpgfb => basis_set_b%pgf_radius, &
225  set_radius_b => basis_set_b%set_radius, &
226  sphi_b => basis_set_b%sphi, &
227  zetb => basis_set_b%zet)
228 
229  nseta = basis_set_a%nset
230  nsetb = basis_set_b%nset
231 
232  IF (inode == 1) last_jatom = 0
233 
234  ! this guarantees minimum image convention
235  ! anything else would not make sense
236  IF (jatom == last_jatom) THEN
237  cycle
238  END IF
239 
240  last_jatom = jatom
241 
242  irow = iatom
243  icol = jatom
244 
245  DO i = 1, 3
246  DO j = 1, 3
247  NULLIFY (mint(i, j)%block)
248  CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
249  row=irow, col=icol, block=mint(i, j)%block, &
250  found=found)
251  cpassert(found)
252  mint(i, j)%block = 0._dp
253  END DO
254  END DO
255 
256  ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
257  ra = pbc(particle_set(iatom)%r(:), cell)
258  rb(:) = ra(:) + rab(:)
259  rac = pbc(rc, ra, cell)
260  rbc = rac + rab
261  dab = norm2(rab)
262 
263  DO iset = 1, nseta
264  ncoa = npgfa(iset)*ncoset(la_max(iset))
265  sgfa = first_sgfa(1, iset)
266  DO jset = 1, nsetb
267  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
268  ncob = npgfb(jset)*ncoset(lb_max(jset))
269  sgfb = first_sgfb(1, jset)
270  ldab = max(ncoa, ncob)
271  lda = ncoset(la_max(iset))*npgfa(iset)
272  ldb = ncoset(lb_max(jset))*npgfb(jset)
273  ALLOCATE (difmab(lda, ldb, m_dim, 3))
274 
275  ! Calculate integral difmab(beta, alpha) = (a|r_beta|db_alpha)
276  ! difmab(beta, alpha) = < a | r_beta | ∂_alpha b >
277  ! difmab(j, idir) = < a | r_j | ∂_idir b >
278  CALL diff_momop_velocity(la_max(iset), npgfa(iset), zeta(:, iset), &
279  rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
280  zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
281  difmab, lambda=lambda, iatom=iatom, jatom=jatom)
282 
283  ! *** Contraction step ***
284 
285  DO idir = 1, 3 ! derivative of AO function
286  DO j = 1, 3 ! position operator r_j
287  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
288  1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
289  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
290  0.0_dp, work(1, 1), SIZE(work, 1))
291 
292  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
293  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
294  work(1, 1), SIZE(work, 1), &
295  1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
296  SIZE(mint(j, idir)%block, 1))
297  END DO !j
298  END DO !idir
299  DEALLOCATE (difmab)
300  END DO !jset
301  END DO !iset
302  END associate
303  END do!iterator
304 
305  CALL neighbor_list_iterator_release(nl_iterator)
306 
307  DO i = 1, 3
308  DO j = 1, 3
309  NULLIFY (mint(i, j)%block)
310  END DO
311  END DO
312 
313  DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
314 
315  CALL timestop(handle)
316  END SUBROUTINE dipole_velocity_deriv
317 
318 ! **************************************************************************************************
319 !> \brief Builds the moments for the derivative of the overlap with respect to nuclear velocities
320 !> \param qs_env ...
321 !> \param moments ...
322 !> \param nmoments ...
323 !> \param ref_point ...
324 !> \param ref_points ...
325 !> \param basis_type ...
326 !> \author Edward Ditler
327 ! **************************************************************************************************
328  SUBROUTINE build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
329 
330  TYPE(qs_environment_type), POINTER :: qs_env
331  TYPE(dbcsr_p_type), DIMENSION(:) :: moments
332  INTEGER, INTENT(IN) :: nmoments
333  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
334  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
335  OPTIONAL :: ref_points
336  CHARACTER(len=*), OPTIONAL :: basis_type
337 
338  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dsdv_moments'
339 
340  INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
341  maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
342  INTEGER, DIMENSION(3) :: image_cell
343  LOGICAL :: found
344  REAL(kind=dp) :: dab
345  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
346  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
347  REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
348  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
349  TYPE(cell_type), POINTER :: cell
350  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
351  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
352  TYPE(neighbor_list_iterator_p_type), &
353  DIMENSION(:), POINTER :: nl_iterator
354  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
355  POINTER :: sab_orb
356  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
357  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
358  TYPE(qs_kind_type), POINTER :: qs_kind
359 
360  IF (nmoments < 1) RETURN
361 
362  CALL timeset(routinen, handle)
363 
364  NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
365 
366  nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
367  cpassert(SIZE(moments) >= nm)
368 
369  NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
370  CALL get_qs_env(qs_env=qs_env, &
371  qs_kind_set=qs_kind_set, &
372  particle_set=particle_set, cell=cell, &
373  sab_orb=sab_orb)
374 
375  nkind = SIZE(qs_kind_set)
376  natom = SIZE(particle_set)
377 
378  ! Allocate work storage
379  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
380  maxco=maxco, maxsgf=maxsgf, &
381  basis_type=basis_type)
382 
383  ALLOCATE (mab(maxco, maxco, nm))
384  mab(:, :, :) = 0.0_dp
385 
386  ALLOCATE (work(maxco, maxsgf))
387  work(:, :) = 0.0_dp
388 
389  ALLOCATE (mint(nm))
390  DO i = 1, nm
391  NULLIFY (mint(i)%block)
392  END DO
393 
394  ALLOCATE (basis_set_list(nkind))
395  DO ikind = 1, nkind
396  qs_kind => qs_kind_set(ikind)
397  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
398  IF (ASSOCIATED(basis_set_a)) THEN
399  basis_set_list(ikind)%gto_basis_set => basis_set_a
400  ELSE
401  NULLIFY (basis_set_list(ikind)%gto_basis_set)
402  END IF
403  END DO
404  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
405  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
406  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
407  iatom=iatom, jatom=jatom, r=rab, cell=image_cell)
408  basis_set_a => basis_set_list(ikind)%gto_basis_set
409  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
410  basis_set_b => basis_set_list(jkind)%gto_basis_set
411  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
412  ! basis ikind
413  associate( &
414  first_sgfa => basis_set_a%first_sgf, &
415  la_max => basis_set_a%lmax, &
416  la_min => basis_set_a%lmin, &
417  npgfa => basis_set_a%npgf, &
418  nsgfa => basis_set_a%nsgf_set, &
419  rpgfa => basis_set_a%pgf_radius, &
420  set_radius_a => basis_set_a%set_radius, &
421  sphi_a => basis_set_a%sphi, &
422  zeta => basis_set_a%zet, &
423  ! basis jkind, &
424  first_sgfb => basis_set_b%first_sgf, &
425  lb_max => basis_set_b%lmax, &
426  lb_min => basis_set_b%lmin, &
427  npgfb => basis_set_b%npgf, &
428  nsgfb => basis_set_b%nsgf_set, &
429  rpgfb => basis_set_b%pgf_radius, &
430  set_radius_b => basis_set_b%set_radius, &
431  sphi_b => basis_set_b%sphi, &
432  zetb => basis_set_b%zet)
433 
434  nseta = basis_set_a%nset
435  nsetb = basis_set_b%nset
436 
437  IF (inode == 1) last_jatom = 0
438 
439  ! this guarantees minimum image convention
440  ! anything else would not make sense
441  IF (jatom == last_jatom) THEN
442  cycle
443  END IF
444 
445  last_jatom = jatom
446 
447  IF (iatom <= jatom) THEN
448  irow = iatom
449  icol = jatom
450  ELSE
451  irow = jatom
452  icol = iatom
453  END IF
454 
455  DO i = 1, nm
456  NULLIFY (mint(i)%block)
457  CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
458  row=irow, col=icol, block=mint(i)%block, found=found)
459  mint(i)%block = 0._dp
460  END DO
461 
462  ! fold atomic position back into unit cell
463  IF (PRESENT(ref_points)) THEN
464  rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
465  ELSE IF (PRESENT(ref_point)) THEN
466  rc(:) = ref_point(:)
467  ELSE
468  rc(:) = 0._dp
469  END IF
470  ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
471  ! by folding around the center, such screwing can be avoided for a proper choice of center.
472  ! we dont use PBC at this point
473 
474  ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
475  ra(:) = particle_set(iatom)%r(:)
476  rb(:) = ra(:) + rab(:)
477  rac = pbc(rc, ra, cell)
478  rbc = rac + rab
479 
480  dab = norm2(rab)
481 
482  DO iset = 1, nseta
483 
484  ncoa = npgfa(iset)*ncoset(la_max(iset))
485  sgfa = first_sgfa(1, iset)
486 
487  DO jset = 1, nsetb
488 
489  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
490 
491  ncob = npgfb(jset)*ncoset(lb_max(jset))
492  sgfb = first_sgfb(1, jset)
493 
494  ! Calculate the primitive integrals
495  CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
496  rpgfa(:, iset), la_min(iset), &
497  lb_max(jset), npgfb(jset), zetb(:, jset), &
498  rpgfb(:, jset), nmoments, rac, rbc, mab)
499 
500  ! Contraction step
501  DO i = 1, nm
502 
503  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
504  1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
505  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
506  0.0_dp, work(1, 1), SIZE(work, 1))
507 
508  IF (iatom <= jatom) THEN
509 
510  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
511  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
512  work(1, 1), SIZE(work, 1), &
513  1.0_dp, mint(i)%block(sgfa, sgfb), &
514  SIZE(mint(i)%block, 1))
515 
516  ELSE
517 
518  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
519  1.0_dp, work(1, 1), SIZE(work, 1), &
520  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
521  1.0_dp, mint(i)%block(sgfb, sgfa), &
522  SIZE(mint(i)%block, 1))
523 
524  END IF
525 
526  END DO
527 
528  END DO
529  END DO
530  END associate
531 
532  END DO ! iterator
533 
534  CALL neighbor_list_iterator_release(nl_iterator)
535 
536  ! Release work storage
537  DEALLOCATE (mab, basis_set_list)
538  DEALLOCATE (work)
539  DO i = 1, nm
540  NULLIFY (mint(i)%block)
541  END DO
542  DEALLOCATE (mint)
543 
544  CALL timestop(handle)
545 
546  END SUBROUTINE build_dsdv_moments
547 
548 ! **************************************************************************************************
549 !> \brief ...
550 !> \param qs_env ...
551 !> \param moments ...
552 !> \param nmoments ...
553 !> \param ref_point ...
554 !> \param ref_points ...
555 !> \param basis_type ...
556 ! **************************************************************************************************
557  SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
558 
559  TYPE(qs_environment_type), POINTER :: qs_env
560  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
561  INTEGER, INTENT(IN) :: nmoments
562  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
563  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
564  OPTIONAL :: ref_points
565  CHARACTER(len=*), OPTIONAL :: basis_type
566 
567  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
568 
569  INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
570  maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
571  LOGICAL :: found
572  REAL(kind=dp) :: dab
573  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
574  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
575  REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
576  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
577  TYPE(cell_type), POINTER :: cell
578  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
579  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
580  TYPE(neighbor_list_iterator_p_type), &
581  DIMENSION(:), POINTER :: nl_iterator
582  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
583  POINTER :: sab_orb
584  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
585  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
586  TYPE(qs_kind_type), POINTER :: qs_kind
587 
588  IF (nmoments < 1) RETURN
589 
590  CALL timeset(routinen, handle)
591 
592  NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
593 
594  nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
595  cpassert(SIZE(moments) >= nm)
596 
597  NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
598  CALL get_qs_env(qs_env=qs_env, &
599  qs_kind_set=qs_kind_set, &
600  particle_set=particle_set, cell=cell, &
601  sab_orb=sab_orb)
602 
603  nkind = SIZE(qs_kind_set)
604  natom = SIZE(particle_set)
605 
606  ! Allocate work storage
607  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
608  maxco=maxco, maxsgf=maxsgf, &
609  basis_type=basis_type)
610 
611  ALLOCATE (mab(maxco, maxco, nm))
612  mab(:, :, :) = 0.0_dp
613 
614  ALLOCATE (work(maxco, maxsgf))
615  work(:, :) = 0.0_dp
616 
617  ALLOCATE (mint(nm))
618  DO i = 1, nm
619  NULLIFY (mint(i)%block)
620  END DO
621 
622  ALLOCATE (basis_set_list(nkind))
623  DO ikind = 1, nkind
624  qs_kind => qs_kind_set(ikind)
625  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
626  IF (ASSOCIATED(basis_set_a)) THEN
627  basis_set_list(ikind)%gto_basis_set => basis_set_a
628  ELSE
629  NULLIFY (basis_set_list(ikind)%gto_basis_set)
630  END IF
631  END DO
632  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
633  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
634  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
635  iatom=iatom, jatom=jatom, r=rab)
636  basis_set_a => basis_set_list(ikind)%gto_basis_set
637  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
638  basis_set_b => basis_set_list(jkind)%gto_basis_set
639  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
640  associate( &
641  ! basis ikind
642  first_sgfa => basis_set_a%first_sgf, &
643  la_max => basis_set_a%lmax, &
644  la_min => basis_set_a%lmin, &
645  npgfa => basis_set_a%npgf, &
646  nsgfa => basis_set_a%nsgf_set, &
647  rpgfa => basis_set_a%pgf_radius, &
648  set_radius_a => basis_set_a%set_radius, &
649  sphi_a => basis_set_a%sphi, &
650  zeta => basis_set_a%zet, &
651  ! basis jkind, &
652  first_sgfb => basis_set_b%first_sgf, &
653  lb_max => basis_set_b%lmax, &
654  lb_min => basis_set_b%lmin, &
655  npgfb => basis_set_b%npgf, &
656  nsgfb => basis_set_b%nsgf_set, &
657  rpgfb => basis_set_b%pgf_radius, &
658  set_radius_b => basis_set_b%set_radius, &
659  sphi_b => basis_set_b%sphi, &
660  zetb => basis_set_b%zet)
661 
662  nseta = basis_set_a%nset
663  nsetb = basis_set_b%nset
664 
665  IF (inode == 1) last_jatom = 0
666 
667  ! this guarantees minimum image convention
668  ! anything else would not make sense
669  IF (jatom == last_jatom) THEN
670  cycle
671  END IF
672 
673  last_jatom = jatom
674 
675  IF (iatom <= jatom) THEN
676  irow = iatom
677  icol = jatom
678  ELSE
679  irow = jatom
680  icol = iatom
681  END IF
682 
683  DO i = 1, nm
684  NULLIFY (mint(i)%block)
685  CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
686  row=irow, col=icol, block=mint(i)%block, found=found)
687  mint(i)%block = 0._dp
688  END DO
689 
690  ! fold atomic position back into unit cell
691  IF (PRESENT(ref_points)) THEN
692  rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
693  ELSE IF (PRESENT(ref_point)) THEN
694  rc(:) = ref_point(:)
695  ELSE
696  rc(:) = 0._dp
697  END IF
698  ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
699  ! by folding around the center, such screwing can be avoided for a proper choice of center.
700  ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
701  rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
702  ! we dont use PBC at this point
703  rab(:) = ra(:) - rb(:)
704  rac(:) = ra(:) - rc(:)
705  rbc(:) = rb(:) - rc(:)
706  dab = norm2(rab)
707 
708  DO iset = 1, nseta
709 
710  ncoa = npgfa(iset)*ncoset(la_max(iset))
711  sgfa = first_sgfa(1, iset)
712 
713  DO jset = 1, nsetb
714 
715  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
716 
717  ncob = npgfb(jset)*ncoset(lb_max(jset))
718  sgfb = first_sgfb(1, jset)
719 
720  ! Calculate the primitive integrals
721  CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
722  rpgfa(:, iset), la_min(iset), &
723  lb_max(jset), npgfb(jset), zetb(:, jset), &
724  rpgfb(:, jset), nmoments, rac, rbc, mab)
725 
726  ! Contraction step
727  DO i = 1, nm
728 
729  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
730  1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
731  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
732  0.0_dp, work(1, 1), SIZE(work, 1))
733 
734  IF (iatom <= jatom) THEN
735 
736  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
737  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
738  work(1, 1), SIZE(work, 1), &
739  1.0_dp, mint(i)%block(sgfa, sgfb), &
740  SIZE(mint(i)%block, 1))
741 
742  ELSE
743 
744  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
745  1.0_dp, work(1, 1), SIZE(work, 1), &
746  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
747  1.0_dp, mint(i)%block(sgfb, sgfa), &
748  SIZE(mint(i)%block, 1))
749 
750  END IF
751 
752  END DO
753 
754  END DO
755  END DO
756  END associate
757 
758  END DO
759  CALL neighbor_list_iterator_release(nl_iterator)
760 
761  ! Release work storage
762  DEALLOCATE (mab, basis_set_list)
763  DEALLOCATE (work)
764  DO i = 1, nm
765  NULLIFY (mint(i)%block)
766  END DO
767  DEALLOCATE (mint)
768 
769  CALL timestop(handle)
770 
771  END SUBROUTINE build_local_moment_matrix
772 
773 ! **************************************************************************************************
774 !> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
775 !> Optionally stores the multipole moments themselves for free.
776 !> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
777 !> Only first derivatives are performed, e. g. x d/dy
778 !> \param qs_env ...
779 !> \param moments_der will contain the derivatives of the multipole moments
780 !> \param nmoments_der order of the moments with derivatives
781 !> \param nmoments order of the multipole moments (no derivatives, same output as
782 !> build_local_moment_matrix, needs moments as arguments to store results)
783 !> \param ref_point ...
784 !> \param moments contains the multipole moments, optionally for free, up to order nmoments
785 !> \note
786 !> Adapted from rRc_xyz_der_ao in qs_operators_ao
787 ! **************************************************************************************************
788  SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
789  ref_point, moments)
790  TYPE(qs_environment_type), POINTER :: qs_env
791  TYPE(dbcsr_p_type), DIMENSION(:, :), &
792  INTENT(INOUT), POINTER :: moments_der
793  INTEGER, INTENT(IN) :: nmoments_der, nmoments
794  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
795  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
796  OPTIONAL, POINTER :: moments
797 
798  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
799 
800  INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
801  jatom, jkind, jpgf, jset, last_jatom, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
802  nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
803  LOGICAL :: found
804  REAL(kind=dp) :: dab
805  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
806  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
807  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
808  REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
809  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
810  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
811  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
812  TYPE(cell_type), POINTER :: cell
813  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
814  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
815  TYPE(neighbor_list_iterator_p_type), &
816  DIMENSION(:), POINTER :: nl_iterator
817  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
818  POINTER :: sab_orb
819  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
820  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
821  TYPE(qs_kind_type), POINTER :: qs_kind
822 
823  nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
824  IF (nmom_build < 1) RETURN
825 
826  CALL timeset(routinen, handle)
827 
828  nders = 1 ! only first order derivatives
829  dimders = ncoset(nders) - 1
830 
831  NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
832  CALL get_qs_env(qs_env=qs_env, &
833  qs_kind_set=qs_kind_set, &
834  particle_set=particle_set, &
835  cell=cell, &
836  sab_orb=sab_orb)
837 
838  nkind = SIZE(qs_kind_set)
839 
840  ! Work storage
841  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
842  maxco=maxco, maxsgf=maxsgf)
843 
844  IF (nmoments > 0) THEN
845  cpassert(PRESENT(moments))
846  nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
847  cpassert(SIZE(moments) == nm)
848  ! storage for integrals
849  ALLOCATE (mab(maxco, maxco, nm))
850  ! blocks
851  mab(:, :, :) = 0.0_dp
852  ALLOCATE (mom_block(nm))
853  DO i = 1, nm
854  NULLIFY (mom_block(i)%block)
855  END DO
856  END IF
857 
858  IF (nmoments_der > 0) THEN
859  m_dim = ncoset(nmoments_der) - 1
860  cpassert(SIZE(moments_der, dim=1) == m_dim)
861  cpassert(SIZE(moments_der, dim=2) == dimders)
862  ! storage for integrals
863  ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
864  difmab(:, :, :, :) = 0.0_dp
865  ! blocks
866  ALLOCATE (mom_block_der(m_dim, dimders))
867  DO i = 1, m_dim
868  DO ider = 1, dimders
869  NULLIFY (mom_block_der(i, ider)%block)
870  END DO
871  END DO
872  END IF
873 
874  ALLOCATE (work(maxco, maxsgf))
875  work(:, :) = 0.0_dp
876 
877  NULLIFY (basis_set_a, basis_set_b, basis_set_list)
878  NULLIFY (qs_kind)
879  ALLOCATE (basis_set_list(nkind))
880  DO ikind = 1, nkind
881  qs_kind => qs_kind_set(ikind)
882  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
883  IF (ASSOCIATED(basis_set_a)) THEN
884  basis_set_list(ikind)%gto_basis_set => basis_set_a
885  ELSE
886  NULLIFY (basis_set_list(ikind)%gto_basis_set)
887  END IF
888  END DO
889 
890  ! Calculate derivatives looping over neighbour list
891  NULLIFY (nl_iterator)
892  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
893  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
894  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
895  iatom=iatom, jatom=jatom, r=rab)
896  basis_set_a => basis_set_list(ikind)%gto_basis_set
897  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
898  basis_set_b => basis_set_list(jkind)%gto_basis_set
899  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
900  associate( &
901  ! basis ikind
902  first_sgfa => basis_set_a%first_sgf, &
903  la_max => basis_set_a%lmax, &
904  la_min => basis_set_a%lmin, &
905  npgfa => basis_set_a%npgf, &
906  nsgfa => basis_set_a%nsgf_set, &
907  rpgfa => basis_set_a%pgf_radius, &
908  set_radius_a => basis_set_a%set_radius, &
909  sphi_a => basis_set_a%sphi, &
910  zeta => basis_set_a%zet, &
911  ! basis jkind, &
912  first_sgfb => basis_set_b%first_sgf, &
913  lb_max => basis_set_b%lmax, &
914  lb_min => basis_set_b%lmin, &
915  npgfb => basis_set_b%npgf, &
916  nsgfb => basis_set_b%nsgf_set, &
917  rpgfb => basis_set_b%pgf_radius, &
918  set_radius_b => basis_set_b%set_radius, &
919  sphi_b => basis_set_b%sphi, &
920  zetb => basis_set_b%zet)
921 
922  nseta = basis_set_a%nset
923  nsetb = basis_set_b%nset
924 
925  ! reference point
926  IF (PRESENT(ref_point)) THEN
927  rc(:) = ref_point(:)
928  ELSE
929  rc(:) = 0._dp
930  END IF
931  ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
932  ! by folding around the center, such screwing can be avoided for a proper choice of center.
933  ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
934  rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
935  ! we dont use PBC at this point
936  rab(:) = ra(:) - rb(:)
937  rac(:) = ra(:) - rc(:)
938  rbc(:) = rb(:) - rc(:)
939  dab = norm2(rab)
940 
941  ! get blocks
942  IF (inode == 1) last_jatom = 0
943 
944  IF (jatom == last_jatom) THEN
945  cycle
946  END IF
947 
948  last_jatom = jatom
949 
950  IF (iatom <= jatom) THEN
951  irow = iatom
952  icol = jatom
953  ELSE
954  irow = jatom
955  icol = iatom
956  END IF
957 
958  IF (nmoments > 0) THEN
959  DO i = 1, nm
960  NULLIFY (mom_block(i)%block)
961  ! get block from pre calculated overlap matrix
962  CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
963  row=irow, col=icol, block=mom_block(i)%block, found=found)
964  cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
965  mom_block(i)%block = 0._dp
966  END DO
967  END IF
968  IF (nmoments_der > 0) THEN
969  DO i = 1, m_dim
970  DO ider = 1, dimders
971  NULLIFY (mom_block_der(i, ider)%block)
972  CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
973  row=irow, col=icol, &
974  block=mom_block_der(i, ider)%block, &
975  found=found)
976  cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
977  mom_block_der(i, ider)%block = 0._dp
978  END DO
979  END DO
980  END IF
981 
982  DO iset = 1, nseta
983 
984  ncoa = npgfa(iset)*ncoset(la_max(iset))
985  sgfa = first_sgfa(1, iset)
986 
987  DO jset = 1, nsetb
988 
989  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
990 
991  ncob = npgfb(jset)*ncoset(lb_max(jset))
992  sgfb = first_sgfb(1, jset)
993 
994  NULLIFY (mab_tmp)
995  ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
996  npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
997 
998  ! Calculate the primitive integrals (need l+1 for derivatives)
999  CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1000  rpgfa(:, iset), la_min(iset), &
1001  lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
1002  rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
1003 
1004  IF (nmoments_der > 0) THEN
1005  CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1006  rpgfa(:, iset), la_min(iset), &
1007  lb_max(jset), npgfb(jset), zetb(:, jset), &
1008  rpgfb(:, jset), lb_min(jset), &
1009  nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp)
1010  END IF
1011 
1012  IF (nmoments > 0) THEN
1013  ! copy subset of mab_tmp (l+1) to mab (l)
1014  mab = 0.0_dp
1015  DO ii = 1, nm
1016  na = 0
1017  nda = 0
1018  DO ipgf = 1, npgfa(iset)
1019  nb = 0
1020  ndb = 0
1021  DO jpgf = 1, npgfb(jset)
1022  DO j = 1, ncoset(lb_max(jset))
1023  DO i = 1, ncoset(la_max(iset))
1024  mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
1025  END DO ! i
1026  END DO ! j
1027  nb = nb + ncoset(lb_max(jset))
1028  ndb = ndb + ncoset(lb_max(jset) + 1)
1029  END DO ! jpgf
1030  na = na + ncoset(la_max(iset))
1031  nda = nda + ncoset(la_max(iset) + 1)
1032  END DO ! ipgf
1033  END DO
1034  ! Contraction step
1035  DO i = 1, nm
1036 
1037  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1038  1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1039  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1040  0.0_dp, work(1, 1), SIZE(work, 1))
1041 
1042  IF (iatom <= jatom) THEN
1043  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1044  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1045  work(1, 1), SIZE(work, 1), &
1046  1.0_dp, mom_block(i)%block(sgfa, sgfb), &
1047  SIZE(mom_block(i)%block, 1))
1048  ELSE
1049  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1050  1.0_dp, work(1, 1), SIZE(work, 1), &
1051  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1052  1.0_dp, mom_block(i)%block(sgfb, sgfa), &
1053  SIZE(mom_block(i)%block, 1))
1054  END IF
1055  END DO
1056  END IF
1057 
1058  IF (nmoments_der > 0) THEN
1059  DO i = 1, m_dim
1060  DO ider = 1, dimders
1061  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1062  1.0_dp, difmab(1, 1, i, ider), SIZE(difmab, 1), &
1063  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1064  0._dp, work(1, 1), SIZE(work, 1))
1065 
1066  IF (iatom <= jatom) THEN
1067  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1068  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1069  work(1, 1), SIZE(work, 1), &
1070  1._dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
1071  SIZE(mom_block_der(i, ider)%block, 1))
1072  ELSE
1073  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1074  -1.0_dp, work(1, 1), SIZE(work, 1), &
1075  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1076  1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
1077  SIZE(mom_block_der(i, ider)%block, 1))
1078  END IF
1079  END DO
1080  END DO
1081  END IF
1082  DEALLOCATE (mab_tmp)
1083  END DO
1084  END DO
1085  END associate
1086  END DO
1087  CALL neighbor_list_iterator_release(nl_iterator)
1088 
1089  ! deallocations
1090  DEALLOCATE (basis_set_list)
1091  DEALLOCATE (work)
1092  IF (nmoments > 0) THEN
1093  DEALLOCATE (mab)
1094  DO i = 1, nm
1095  NULLIFY (mom_block(i)%block)
1096  END DO
1097  DEALLOCATE (mom_block)
1098  END IF
1099  IF (nmoments_der > 0) THEN
1100  DEALLOCATE (difmab)
1101  DO i = 1, m_dim
1102  DO ider = 1, dimders
1103  NULLIFY (mom_block_der(i, ider)%block)
1104  END DO
1105  END DO
1106  DEALLOCATE (mom_block_der)
1107  END IF
1108 
1109  CALL timestop(handle)
1110 
1111  END SUBROUTINE build_local_moments_der_matrix
1112 
1113 ! **************************************************************************************************
1114 !> \brief ...
1115 !> \param qs_env ...
1116 !> \param magmom ...
1117 !> \param nmoments ...
1118 !> \param ref_point ...
1119 !> \param ref_points ...
1120 !> \param basis_type ...
1121 ! **************************************************************************************************
1122  SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
1123 
1124  TYPE(qs_environment_type), POINTER :: qs_env
1125  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
1126  INTEGER, INTENT(IN) :: nmoments
1127  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
1128  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1129  OPTIONAL :: ref_points
1130  CHARACTER(len=*), OPTIONAL :: basis_type
1131 
1132  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
1133 
1134  INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
1135  maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
1136  LOGICAL :: found
1137  REAL(kind=dp) :: dab
1138  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1139  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
1140  REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
1141  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
1142  TYPE(cell_type), POINTER :: cell
1143  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1144  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1145  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1146  TYPE(neighbor_list_iterator_p_type), &
1147  DIMENSION(:), POINTER :: nl_iterator
1148  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1149  POINTER :: sab_orb
1150  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1151  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1152  TYPE(qs_kind_type), POINTER :: qs_kind
1153 
1154  IF (nmoments < 1) RETURN
1155 
1156  CALL timeset(routinen, handle)
1157 
1158  NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
1159  NULLIFY (matrix_s)
1160 
1161  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1162 
1163  ! magnetic dipoles/angular moments only
1164  nm = 3
1165 
1166  NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1167  CALL get_qs_env(qs_env=qs_env, &
1168  qs_kind_set=qs_kind_set, &
1169  particle_set=particle_set, cell=cell, &
1170  sab_orb=sab_orb)
1171 
1172  nkind = SIZE(qs_kind_set)
1173  natom = SIZE(particle_set)
1174 
1175  ! Allocate work storage
1176  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1177  maxco=maxco, maxsgf=maxsgf)
1178 
1179  ALLOCATE (mab(maxco, maxco, nm))
1180  mab(:, :, :) = 0.0_dp
1181 
1182  ALLOCATE (work(maxco, maxsgf))
1183  work(:, :) = 0.0_dp
1184 
1185  ALLOCATE (mint(nm))
1186  DO i = 1, nm
1187  NULLIFY (mint(i)%block)
1188  END DO
1189 
1190  ALLOCATE (basis_set_list(nkind))
1191  DO ikind = 1, nkind
1192  qs_kind => qs_kind_set(ikind)
1193  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1194  IF (ASSOCIATED(basis_set_a)) THEN
1195  basis_set_list(ikind)%gto_basis_set => basis_set_a
1196  ELSE
1197  NULLIFY (basis_set_list(ikind)%gto_basis_set)
1198  END IF
1199  END DO
1200  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1201  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1202  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1203  iatom=iatom, jatom=jatom, r=rab)
1204  basis_set_a => basis_set_list(ikind)%gto_basis_set
1205  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1206  basis_set_b => basis_set_list(jkind)%gto_basis_set
1207  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1208  associate( &
1209  ! basis ikind
1210  first_sgfa => basis_set_a%first_sgf, &
1211  la_max => basis_set_a%lmax, &
1212  la_min => basis_set_a%lmin, &
1213  npgfa => basis_set_a%npgf, &
1214  nsgfa => basis_set_a%nsgf_set, &
1215  rpgfa => basis_set_a%pgf_radius, &
1216  set_radius_a => basis_set_a%set_radius, &
1217  sphi_a => basis_set_a%sphi, &
1218  zeta => basis_set_a%zet, &
1219  ! basis jkind, &
1220  first_sgfb => basis_set_b%first_sgf, &
1221  lb_max => basis_set_b%lmax, &
1222  lb_min => basis_set_b%lmin, &
1223  npgfb => basis_set_b%npgf, &
1224  nsgfb => basis_set_b%nsgf_set, &
1225  rpgfb => basis_set_b%pgf_radius, &
1226  set_radius_b => basis_set_b%set_radius, &
1227  sphi_b => basis_set_b%sphi, &
1228  zetb => basis_set_b%zet)
1229 
1230  nseta = basis_set_a%nset
1231  nsetb = basis_set_b%nset
1232 
1233  IF (iatom <= jatom) THEN
1234  irow = iatom
1235  icol = jatom
1236  ELSE
1237  irow = jatom
1238  icol = iatom
1239  END IF
1240 
1241  DO i = 1, nm
1242  NULLIFY (mint(i)%block)
1243  CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1244  row=irow, col=icol, block=mint(i)%block, found=found)
1245  mint(i)%block = 0._dp
1246  cpassert(ASSOCIATED(mint(i)%block))
1247  END DO
1248 
1249  ! fold atomic position back into unit cell
1250  IF (PRESENT(ref_points)) THEN
1251  rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1252  ELSE IF (PRESENT(ref_point)) THEN
1253  rc(:) = ref_point(:)
1254  ELSE
1255  rc(:) = 0._dp
1256  END IF
1257  ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1258  ! by folding around the center, such screwing can be avoided for a proper choice of center.
1259  ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1260  rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1261  ! we dont use PBC at this point
1262  rab(:) = ra(:) - rb(:)
1263  rac(:) = ra(:) - rc(:)
1264  rbc(:) = rb(:) - rc(:)
1265  dab = norm2(rab)
1266 
1267  DO iset = 1, nseta
1268 
1269  ncoa = npgfa(iset)*ncoset(la_max(iset))
1270  sgfa = first_sgfa(1, iset)
1271 
1272  DO jset = 1, nsetb
1273 
1274  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1275 
1276  ncob = npgfb(jset)*ncoset(lb_max(jset))
1277  sgfb = first_sgfb(1, jset)
1278 
1279  ! Calculate the primitive integrals
1280  CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1281  rpgfa(:, iset), la_min(iset), &
1282  lb_max(jset), npgfb(jset), zetb(:, jset), &
1283  rpgfb(:, jset), rac, rbc, mab)
1284 
1285  ! Contraction step
1286  DO i = 1, nm
1287  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1288  1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1289  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1290  0.0_dp, work(1, 1), SIZE(work, 1))
1291 
1292  IF (iatom <= jatom) THEN
1293  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1294  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1295  work(1, 1), SIZE(work, 1), &
1296  1.0_dp, mint(i)%block(sgfa, sgfb), &
1297  SIZE(mint(i)%block, 1))
1298  ELSE
1299  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1300  -1.0_dp, work(1, 1), SIZE(work, 1), &
1301  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1302  1.0_dp, mint(i)%block(sgfb, sgfa), &
1303  SIZE(mint(i)%block, 1))
1304  END IF
1305 
1306  END DO
1307 
1308  END DO
1309  END DO
1310  END associate
1311  END DO
1312  CALL neighbor_list_iterator_release(nl_iterator)
1313 
1314  ! Release work storage
1315  DEALLOCATE (mab, basis_set_list)
1316  DEALLOCATE (work)
1317  DO i = 1, nm
1318  NULLIFY (mint(i)%block)
1319  END DO
1320  DEALLOCATE (mint)
1321 
1322  CALL timestop(handle)
1323 
1324  END SUBROUTINE build_local_magmom_matrix
1325 
1326 ! **************************************************************************************************
1327 !> \brief ...
1328 !> \param qs_env ...
1329 !> \param cosmat ...
1330 !> \param sinmat ...
1331 !> \param kvec ...
1332 !> \param sab_orb_external ...
1333 !> \param basis_type ...
1334 ! **************************************************************************************************
1335  SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1336 
1337  TYPE(qs_environment_type), POINTER :: qs_env
1338  TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1339  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1340  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1341  OPTIONAL, POINTER :: sab_orb_external
1342  CHARACTER(len=*), OPTIONAL :: basis_type
1343 
1344  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1345 
1346  INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, ldsa, &
1347  ldsb, ldwork, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
1348  LOGICAL :: found
1349  REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1350  REAL(kind=dp) :: dab
1351  REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1352  TYPE(cell_type), POINTER :: cell
1353  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1354  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1355  TYPE(neighbor_list_iterator_p_type), &
1356  DIMENSION(:), POINTER :: nl_iterator
1357  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1358  POINTER :: sab_orb
1359  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1360  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1361  TYPE(qs_kind_type), POINTER :: qs_kind
1362 
1363  CALL timeset(routinen, handle)
1364 
1365  NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1366  CALL get_qs_env(qs_env=qs_env, &
1367  qs_kind_set=qs_kind_set, &
1368  particle_set=particle_set, cell=cell, &
1369  sab_orb=sab_orb)
1370 
1371  IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
1372 
1373  CALL dbcsr_set(sinmat, 0.0_dp)
1374  CALL dbcsr_set(cosmat, 0.0_dp)
1375 
1376  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
1377  ldab = ldwork
1378  ALLOCATE (cosab(ldab, ldab))
1379  ALLOCATE (sinab(ldab, ldab))
1380  ALLOCATE (work(ldwork, ldwork))
1381 
1382  nkind = SIZE(qs_kind_set)
1383  natom = SIZE(particle_set)
1384 
1385  ALLOCATE (basis_set_list(nkind))
1386  DO ikind = 1, nkind
1387  qs_kind => qs_kind_set(ikind)
1388  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1389  IF (ASSOCIATED(basis_set_a)) THEN
1390  basis_set_list(ikind)%gto_basis_set => basis_set_a
1391  ELSE
1392  NULLIFY (basis_set_list(ikind)%gto_basis_set)
1393  END IF
1394  END DO
1395  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1396  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1397  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1398  iatom=iatom, jatom=jatom, r=rab)
1399  basis_set_a => basis_set_list(ikind)%gto_basis_set
1400  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1401  basis_set_b => basis_set_list(jkind)%gto_basis_set
1402  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1403  associate( &
1404  ! basis ikind
1405  first_sgfa => basis_set_a%first_sgf, &
1406  la_max => basis_set_a%lmax, &
1407  la_min => basis_set_a%lmin, &
1408  npgfa => basis_set_a%npgf, &
1409  nsgfa => basis_set_a%nsgf_set, &
1410  rpgfa => basis_set_a%pgf_radius, &
1411  set_radius_a => basis_set_a%set_radius, &
1412  sphi_a => basis_set_a%sphi, &
1413  zeta => basis_set_a%zet, &
1414  ! basis jkind, &
1415  first_sgfb => basis_set_b%first_sgf, &
1416  lb_max => basis_set_b%lmax, &
1417  lb_min => basis_set_b%lmin, &
1418  npgfb => basis_set_b%npgf, &
1419  nsgfb => basis_set_b%nsgf_set, &
1420  rpgfb => basis_set_b%pgf_radius, &
1421  set_radius_b => basis_set_b%set_radius, &
1422  sphi_b => basis_set_b%sphi, &
1423  zetb => basis_set_b%zet)
1424 
1425  nseta = basis_set_a%nset
1426  nsetb = basis_set_b%nset
1427 
1428  ldsa = SIZE(sphi_a, 1)
1429  ldsb = SIZE(sphi_b, 1)
1430 
1431  IF (iatom <= jatom) THEN
1432  irow = iatom
1433  icol = jatom
1434  ELSE
1435  irow = jatom
1436  icol = iatom
1437  END IF
1438 
1439  NULLIFY (cblock)
1440  CALL dbcsr_get_block_p(matrix=cosmat, &
1441  row=irow, col=icol, block=cblock, found=found)
1442  NULLIFY (sblock)
1443  CALL dbcsr_get_block_p(matrix=sinmat, &
1444  row=irow, col=icol, block=sblock, found=found)
1445  IF (ASSOCIATED(cblock) .AND. .NOT. ASSOCIATED(sblock) .OR. &
1446  .NOT. ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1447  cpabort("")
1448  END IF
1449 
1450  IF (ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1451 
1452  ra(:) = pbc(particle_set(iatom)%r(:), cell)
1453  rb(:) = ra + rab
1454  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1455 
1456  DO iset = 1, nseta
1457 
1458  ncoa = npgfa(iset)*ncoset(la_max(iset))
1459  sgfa = first_sgfa(1, iset)
1460 
1461  DO jset = 1, nsetb
1462 
1463  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1464 
1465  ncob = npgfb(jset)*ncoset(lb_max(jset))
1466  sgfb = first_sgfb(1, jset)
1467 
1468  ! Calculate the primitive integrals
1469  CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1470  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1471  ra, rb, kvec, cosab, sinab)
1472  CALL contract_cossin(cblock, sblock, &
1473  iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1474  jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1475  cosab, sinab, ldab, work, ldwork)
1476 
1477  END DO
1478  END DO
1479 
1480  END IF
1481  END associate
1482  END DO
1483  CALL neighbor_list_iterator_release(nl_iterator)
1484 
1485  DEALLOCATE (cosab)
1486  DEALLOCATE (sinab)
1487  DEALLOCATE (work)
1488  DEALLOCATE (basis_set_list)
1489 
1490  CALL timestop(handle)
1491 
1492  END SUBROUTINE build_berry_moment_matrix
1493 
1494 ! **************************************************************************************************
1495 !> \brief ...
1496 !> \param qs_env ...
1497 !> \param cosmat ...
1498 !> \param sinmat ...
1499 !> \param kvec ...
1500 !> \param basis_type ...
1501 ! **************************************************************************************************
1502  SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1503 
1504  TYPE(qs_environment_type), POINTER :: qs_env
1505  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1506  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1507  CHARACTER(len=*), OPTIONAL :: basis_type
1508 
1509  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1510 
1511  INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1512  ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1513  INTEGER, DIMENSION(3) :: icell
1514  INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1515  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1516  LOGICAL :: found, use_cell_mapping
1517  REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1518  REAL(kind=dp) :: dab
1519  REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1520  TYPE(cell_type), POINTER :: cell
1521  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1522  TYPE(dft_control_type), POINTER :: dft_control
1523  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1524  TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1525  TYPE(kpoint_type), POINTER :: kpoints
1526  TYPE(neighbor_list_iterator_p_type), &
1527  DIMENSION(:), POINTER :: nl_iterator
1528  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1529  POINTER :: sab_orb
1530  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1531  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1532  TYPE(qs_kind_type), POINTER :: qs_kind
1533  TYPE(qs_ks_env_type), POINTER :: ks_env
1534 
1535  CALL timeset(routinen, handle)
1536 
1537  CALL get_qs_env(qs_env, &
1538  ks_env=ks_env, &
1539  dft_control=dft_control)
1540  nimg = dft_control%nimages
1541  IF (nimg > 1) THEN
1542  CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1543  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1544  use_cell_mapping = .true.
1545  ELSE
1546  use_cell_mapping = .false.
1547  END IF
1548 
1549  CALL get_qs_env(qs_env=qs_env, &
1550  qs_kind_set=qs_kind_set, &
1551  particle_set=particle_set, cell=cell, &
1552  sab_orb=sab_orb)
1553 
1554  nkind = SIZE(qs_kind_set)
1555  natom = SIZE(particle_set)
1556  ALLOCATE (basis_set_list(nkind))
1557  DO ikind = 1, nkind
1558  qs_kind => qs_kind_set(ikind)
1559  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1560  IF (ASSOCIATED(basis_set)) THEN
1561  basis_set_list(ikind)%gto_basis_set => basis_set
1562  ELSE
1563  NULLIFY (basis_set_list(ikind)%gto_basis_set)
1564  END IF
1565  END DO
1566 
1567  ALLOCATE (row_blk_sizes(natom))
1568  CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1569  basis=basis_set_list)
1570  CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1571  ! (re)allocate matrix sets
1572  CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1573  CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1574  DO i = 1, nimg
1575  ! sin
1576  ALLOCATE (sinmat(1, i)%matrix)
1577  CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1578  name="SINMAT", &
1579  dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1580  row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1581  nze=0)
1582  CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1583  CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1584  ! cos
1585  ALLOCATE (cosmat(1, i)%matrix)
1586  CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1587  name="COSMAT", &
1588  dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1589  row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1590  nze=0)
1591  CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1592  CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1593  END DO
1594 
1595  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1596  ldab = ldwork
1597  ALLOCATE (cosab(ldab, ldab))
1598  ALLOCATE (sinab(ldab, ldab))
1599  ALLOCATE (work(ldwork, ldwork))
1600 
1601  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1602  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1603  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1604  iatom=iatom, jatom=jatom, r=rab, cell=icell)
1605  basis_set_a => basis_set_list(ikind)%gto_basis_set
1606  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1607  basis_set_b => basis_set_list(jkind)%gto_basis_set
1608  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1609  associate( &
1610  ! basis ikind
1611  first_sgfa => basis_set_a%first_sgf, &
1612  la_max => basis_set_a%lmax, &
1613  la_min => basis_set_a%lmin, &
1614  npgfa => basis_set_a%npgf, &
1615  nsgfa => basis_set_a%nsgf_set, &
1616  rpgfa => basis_set_a%pgf_radius, &
1617  set_radius_a => basis_set_a%set_radius, &
1618  sphi_a => basis_set_a%sphi, &
1619  zeta => basis_set_a%zet, &
1620  ! basis jkind, &
1621  first_sgfb => basis_set_b%first_sgf, &
1622  lb_max => basis_set_b%lmax, &
1623  lb_min => basis_set_b%lmin, &
1624  npgfb => basis_set_b%npgf, &
1625  nsgfb => basis_set_b%nsgf_set, &
1626  rpgfb => basis_set_b%pgf_radius, &
1627  set_radius_b => basis_set_b%set_radius, &
1628  sphi_b => basis_set_b%sphi, &
1629  zetb => basis_set_b%zet)
1630 
1631  nseta = basis_set_a%nset
1632  nsetb = basis_set_b%nset
1633 
1634  ldsa = SIZE(sphi_a, 1)
1635  ldsb = SIZE(sphi_b, 1)
1636 
1637  IF (iatom <= jatom) THEN
1638  irow = iatom
1639  icol = jatom
1640  ELSE
1641  irow = jatom
1642  icol = iatom
1643  END IF
1644 
1645  IF (use_cell_mapping) THEN
1646  ic = cell_to_index(icell(1), icell(2), icell(3))
1647  cpassert(ic > 0)
1648  ELSE
1649  ic = 1
1650  END IF
1651 
1652  NULLIFY (sblock)
1653  CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1654  row=irow, col=icol, block=sblock, found=found)
1655  cpassert(found)
1656  NULLIFY (cblock)
1657  CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1658  row=irow, col=icol, block=cblock, found=found)
1659  cpassert(found)
1660 
1661  ra(:) = pbc(particle_set(iatom)%r(:), cell)
1662  rb(:) = ra + rab
1663  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1664 
1665  DO iset = 1, nseta
1666 
1667  ncoa = npgfa(iset)*ncoset(la_max(iset))
1668  sgfa = first_sgfa(1, iset)
1669 
1670  DO jset = 1, nsetb
1671 
1672  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1673 
1674  ncob = npgfb(jset)*ncoset(lb_max(jset))
1675  sgfb = first_sgfb(1, jset)
1676 
1677  ! Calculate the primitive integrals
1678  CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1679  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1680  ra, rb, kvec, cosab, sinab)
1681  CALL contract_cossin(cblock, sblock, &
1682  iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1683  jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1684  cosab, sinab, ldab, work, ldwork)
1685 
1686  END DO
1687  END DO
1688  END associate
1689  END DO
1690  CALL neighbor_list_iterator_release(nl_iterator)
1691 
1692  DEALLOCATE (cosab)
1693  DEALLOCATE (sinab)
1694  DEALLOCATE (work)
1695  DEALLOCATE (basis_set_list)
1696  DEALLOCATE (row_blk_sizes)
1697 
1698  CALL timestop(handle)
1699 
1700  END SUBROUTINE build_berry_kpoint_matrix
1701 
1702 ! **************************************************************************************************
1703 !> \brief ...
1704 !> \param qs_env ...
1705 !> \param magnetic ...
1706 !> \param nmoments ...
1707 !> \param reference ...
1708 !> \param ref_point ...
1709 !> \param unit_number ...
1710 ! **************************************************************************************************
1711  SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1712 
1713  TYPE(qs_environment_type), POINTER :: qs_env
1714  LOGICAL, INTENT(IN) :: magnetic
1715  INTEGER, INTENT(IN) :: nmoments, reference
1716  REAL(dp), DIMENSION(:), POINTER :: ref_point
1717  INTEGER, INTENT(IN) :: unit_number
1718 
1719  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_berry_phase'
1720 
1721  CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1722  CHARACTER(LEN=default_string_length) :: description
1723  COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1724  zij(3, 3), zijk(3, 3, 3), &
1725  zijkl(3, 3, 3, 3), zphase(3), zz
1726  INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1727  iy, iz, j, k, l, nao, nm, nmo, nmom, &
1728  nmotot, tmp_dim
1729  LOGICAL :: floating, ghost, uniform
1730  REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1731  REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1732  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1733  REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1734  TYPE(atomic_kind_type), POINTER :: atomic_kind
1735  TYPE(cell_type), POINTER :: cell
1736  TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1737  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1738  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1739  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1740  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1741  TYPE(cp_fm_type), POINTER :: mo_coeff
1742  TYPE(cp_result_type), POINTER :: results
1743  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1744  TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1745  TYPE(dft_control_type), POINTER :: dft_control
1746  TYPE(distribution_1d_type), POINTER :: local_particles
1747  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1748  TYPE(mp_para_env_type), POINTER :: para_env
1749  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1750  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1751  TYPE(qs_rho_type), POINTER :: rho
1752  TYPE(rt_prop_type), POINTER :: rtp
1753 
1754  cpassert(ASSOCIATED(qs_env))
1755 
1756  IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1757  IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1758  RETURN
1759  END IF
1760 
1761  CALL timeset(routinen, handle)
1762 
1763  ! restrict maximum moment available
1764  nmom = min(nmoments, 2)
1765 
1766  nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1767  ! rmom(:,1)=electronic
1768  ! rmom(:,2)=nuclear
1769  ! rmom(:,1)=total
1770  ALLOCATE (rmom(nm + 1, 3))
1771  ALLOCATE (rlab(nm + 1))
1772  rmom = 0.0_dp
1773  rlab = ""
1774  IF (magnetic) THEN
1775  nm = 3
1776  ALLOCATE (mmom(nm))
1777  mmom = 0._dp
1778  END IF
1779 
1780  NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1781  local_particles, matrix_s, mos, rho_ao)
1782 
1783  CALL get_qs_env(qs_env, &
1784  dft_control=dft_control, &
1785  rho=rho, &
1786  cell=cell, &
1787  results=results, &
1788  particle_set=particle_set, &
1789  qs_kind_set=qs_kind_set, &
1790  para_env=para_env, &
1791  local_particles=local_particles, &
1792  matrix_s=matrix_s, &
1793  mos=mos)
1794 
1795  CALL qs_rho_get(rho, rho_ao=rho_ao)
1796 
1797  NULLIFY (cosmat, sinmat)
1798  ALLOCATE (cosmat, sinmat)
1799  CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1800  CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1801  CALL dbcsr_set(cosmat, 0.0_dp)
1802  CALL dbcsr_set(sinmat, 0.0_dp)
1803 
1804  ALLOCATE (op_fm_set(2, dft_control%nspins))
1805  ALLOCATE (opvec(dft_control%nspins))
1806  ALLOCATE (eigrmat(dft_control%nspins))
1807  nmotot = 0
1808  DO ispin = 1, dft_control%nspins
1809  NULLIFY (tmp_fm_struct, mo_coeff)
1810  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1811  nmotot = nmotot + nmo
1812  CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1813  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1814  ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1815  DO i = 1, SIZE(op_fm_set, 1)
1816  CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1817  END DO
1818  CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1819  CALL cp_fm_struct_release(tmp_fm_struct)
1820  END DO
1821 
1822  ! occupation
1823  DO ispin = 1, dft_control%nspins
1824  CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1825  IF (.NOT. uniform) THEN
1826  cpwarn("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1827  END IF
1828  END DO
1829 
1830  ! reference point
1831  CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1832  rcc = pbc(rcc, cell)
1833 
1834  ! label
1835  DO l = 1, nm
1836  ix = indco(1, l + 1)
1837  iy = indco(2, l + 1)
1838  iz = indco(3, l + 1)
1839  CALL set_label(rlab(l + 1), ix, iy, iz)
1840  END DO
1841 
1842  ! nuclear contribution
1843  DO ia = 1, SIZE(particle_set)
1844  atomic_kind => particle_set(ia)%atomic_kind
1845  CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1846  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1847  IF (.NOT. ghost .AND. .NOT. floating) THEN
1848  rmom(1, 2) = rmom(1, 2) - charge
1849  END IF
1850  END DO
1851  ria = twopi*matmul(cell%h_inv, rcc)
1852  zphase = cmplx(cos(ria), sin(ria), dp)**rmom(1, 2)
1853 
1854  zi = 0._dp
1855  zij = 0._dp
1856  zijk = 0._dp
1857  zijkl = 0._dp
1858 
1859  DO l = 1, nmom
1860  SELECT CASE (l)
1861  CASE (1)
1862  ! Dipole
1863  zi(:) = cmplx(1._dp, 0._dp, dp)
1864  DO ia = 1, SIZE(particle_set)
1865  atomic_kind => particle_set(ia)%atomic_kind
1866  CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1867  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1868  IF (.NOT. ghost .AND. .NOT. floating) THEN
1869  ria = particle_set(ia)%r
1870  ria = pbc(ria, cell)
1871  DO i = 1, 3
1872  kvec(:) = twopi*cell%h_inv(i, :)
1873  dd = sum(kvec(:)*ria(:))
1874  zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1875  zi(i) = zi(i)*zdeta
1876  END DO
1877  END IF
1878  END DO
1879  zi = zi*zphase
1880  ci = aimag(log(zi))/twopi
1881  qq = aimag(log(zi))
1882  rmom(2:4, 2) = matmul(cell%hmat, ci)
1883  CASE (2)
1884  ! Quadrupole
1885  cpabort("Berry phase moments bigger than 1 not implemented")
1886  zij(:, :) = cmplx(1._dp, 0._dp, dp)
1887  DO ia = 1, SIZE(particle_set)
1888  atomic_kind => particle_set(ia)%atomic_kind
1889  CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1890  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1891  ria = particle_set(ia)%r
1892  ria = pbc(ria, cell)
1893  DO i = 1, 3
1894  DO j = i, 3
1895  kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1896  dd = sum(kvec(:)*ria(:))
1897  zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1898  zij(i, j) = zij(i, j)*zdeta
1899  zij(j, i) = zij(i, j)
1900  END DO
1901  END DO
1902  END DO
1903  DO i = 1, 3
1904  DO j = 1, 3
1905  zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1906  zz = zij(i, j)/zi(i)/zi(j)
1907  cij(i, j) = aimag(log(zz))/twopi
1908  END DO
1909  END DO
1910  cij = 0.5_dp*cij/twopi/twopi
1911  cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1912  DO k = 4, 9
1913  ix = indco(1, k + 1)
1914  iy = indco(2, k + 1)
1915  iz = indco(3, k + 1)
1916  IF (ix == 0) THEN
1917  rmom(k + 1, 2) = cij(iy, iz)
1918  ELSE IF (iy == 0) THEN
1919  rmom(k + 1, 2) = cij(ix, iz)
1920  ELSE IF (iz == 0) THEN
1921  rmom(k + 1, 2) = cij(ix, iy)
1922  END IF
1923  END DO
1924  CASE (3)
1925  ! Octapole
1926  cpabort("Berry phase moments bigger than 2 not implemented")
1927  CASE (4)
1928  ! Hexadecapole
1929  cpabort("Berry phase moments bigger than 3 not implemented")
1930  CASE DEFAULT
1931  cpabort("Berry phase moments bigger than 4 not implemented")
1932  END SELECT
1933  END DO
1934 
1935  ! electronic contribution
1936 
1937  ria = twopi*real(nmotot, dp)*occ*matmul(cell%h_inv, rcc)
1938  xphase = cmplx(cos(ria), sin(ria), dp)
1939 
1940  ! charge
1941  trace = 0.0_dp
1942  DO ispin = 1, dft_control%nspins
1943  CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1944  rmom(1, 1) = rmom(1, 1) + trace
1945  END DO
1946 
1947  zi = 0._dp
1948  zij = 0._dp
1949  zijk = 0._dp
1950  zijkl = 0._dp
1951 
1952  DO l = 1, nmom
1953  SELECT CASE (l)
1954  CASE (1)
1955  ! Dipole
1956  DO i = 1, 3
1957  kvec(:) = twopi*cell%h_inv(i, :)
1958  CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1959  IF (qs_env%run_rtp) THEN
1960  CALL get_qs_env(qs_env, rtp=rtp)
1961  CALL get_rtp(rtp, mos_new=mos_new)
1962  CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1963  ELSE
1964  CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1965  END IF
1966  zdet = cmplx(1._dp, 0._dp, dp)
1967  DO ispin = 1, dft_control%nspins
1968  CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1969  DO idim = 1, tmp_dim
1970  eigrmat(ispin)%local_data(:, idim) = &
1971  cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1972  -op_fm_set(2, ispin)%local_data(:, idim), dp)
1973  END DO
1974  ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1975  CALL cp_cfm_det(eigrmat(ispin), zdeta)
1976  zdet = zdet*zdeta
1977  IF (dft_control%nspins == 1) THEN
1978  zdet = zdet*zdeta
1979  END IF
1980  END DO
1981  zi(i) = zdet
1982  END DO
1983  zi = zi*xphase
1984  CASE (2)
1985  ! Quadrupole
1986  cpabort("Berry phase moments bigger than 1 not implemented")
1987  DO i = 1, 3
1988  DO j = i, 3
1989  kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1990  CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1991  IF (qs_env%run_rtp) THEN
1992  CALL get_qs_env(qs_env, rtp=rtp)
1993  CALL get_rtp(rtp, mos_new=mos_new)
1994  CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1995  ELSE
1996  CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1997  END IF
1998  zdet = cmplx(1._dp, 0._dp, dp)
1999  DO ispin = 1, dft_control%nspins
2000  CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
2001  DO idim = 1, tmp_dim
2002  eigrmat(ispin)%local_data(:, idim) = &
2003  cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
2004  -op_fm_set(2, ispin)%local_data(:, idim), dp)
2005  END DO
2006  ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
2007  CALL cp_cfm_det(eigrmat(ispin), zdeta)
2008  zdet = zdet*zdeta
2009  IF (dft_control%nspins == 1) THEN
2010  zdet = zdet*zdeta
2011  END IF
2012  END DO
2013  zij(i, j) = zdet*xphase(i)*xphase(j)
2014  zij(j, i) = zdet*xphase(i)*xphase(j)
2015  END DO
2016  END DO
2017  CASE (3)
2018  ! Octapole
2019  cpabort("Berry phase moments bigger than 2 not implemented")
2020  CASE (4)
2021  ! Hexadecapole
2022  cpabort("Berry phase moments bigger than 3 not implemented")
2023  CASE DEFAULT
2024  cpabort("Berry phase moments bigger than 4 not implemented")
2025  END SELECT
2026  END DO
2027  DO l = 1, nmom
2028  SELECT CASE (l)
2029  CASE (1)
2030  ! Dipole (apply periodic (2 Pi) boundary conditions)
2031  ci = aimag(log(zi))
2032  DO i = 1, 3
2033  IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
2034  IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
2035  END DO
2036  rmom(2:4, 1) = matmul(cell%hmat, ci)/twopi
2037  CASE (2)
2038  ! Quadrupole
2039  cpabort("Berry phase moments bigger than 1 not implemented")
2040  DO i = 1, 3
2041  DO j = 1, 3
2042  zz = zij(i, j)/zi(i)/zi(j)
2043  cij(i, j) = aimag(log(zz))/twopi
2044  END DO
2045  END DO
2046  cij = 0.5_dp*cij/twopi/twopi
2047  cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
2048  DO k = 4, 9
2049  ix = indco(1, k + 1)
2050  iy = indco(2, k + 1)
2051  iz = indco(3, k + 1)
2052  IF (ix == 0) THEN
2053  rmom(k + 1, 1) = cij(iy, iz)
2054  ELSE IF (iy == 0) THEN
2055  rmom(k + 1, 1) = cij(ix, iz)
2056  ELSE IF (iz == 0) THEN
2057  rmom(k + 1, 1) = cij(ix, iy)
2058  END IF
2059  END DO
2060  CASE (3)
2061  ! Octapole
2062  cpabort("Berry phase moments bigger than 2 not implemented")
2063  CASE (4)
2064  ! Hexadecapole
2065  cpabort("Berry phase moments bigger than 3 not implemented")
2066  CASE DEFAULT
2067  cpabort("Berry phase moments bigger than 4 not implemented")
2068  END SELECT
2069  END DO
2070 
2071  rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2072  description = "[DIPOLE]"
2073  CALL cp_results_erase(results=results, description=description)
2074  CALL put_results(results=results, description=description, &
2075  values=rmom(2:4, 3))
2076  IF (magnetic) THEN
2077  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., mmom=mmom)
2078  ELSE
2079  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true.)
2080  END IF
2081 
2082  DEALLOCATE (rmom)
2083  DEALLOCATE (rlab)
2084  IF (magnetic) THEN
2085  DEALLOCATE (mmom)
2086  END IF
2087 
2088  CALL dbcsr_deallocate_matrix(cosmat)
2089  CALL dbcsr_deallocate_matrix(sinmat)
2090 
2091  CALL cp_fm_release(opvec)
2092  CALL cp_fm_release(op_fm_set)
2093  DO ispin = 1, dft_control%nspins
2094  CALL cp_cfm_release(eigrmat(ispin))
2095  END DO
2096  DEALLOCATE (eigrmat)
2097 
2098  CALL timestop(handle)
2099 
2100  END SUBROUTINE qs_moment_berry_phase
2101 
2102 ! **************************************************************************************************
2103 !> \brief ...
2104 !> \param cosmat ...
2105 !> \param sinmat ...
2106 !> \param mos ...
2107 !> \param op_fm_set ...
2108 !> \param opvec ...
2109 ! **************************************************************************************************
2110  SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
2111 
2112  TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2113  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2114  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2115  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: opvec
2116 
2117  INTEGER :: i, nao, nmo
2118  TYPE(cp_fm_type), POINTER :: mo_coeff
2119 
2120  DO i = 1, SIZE(op_fm_set, 2) ! spin
2121  CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2122  CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
2123  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2124  op_fm_set(1, i))
2125  CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
2126  CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2127  op_fm_set(2, i))
2128  END DO
2129 
2130  END SUBROUTINE op_orbbas
2131 
2132 ! **************************************************************************************************
2133 !> \brief ...
2134 !> \param cosmat ...
2135 !> \param sinmat ...
2136 !> \param mos ...
2137 !> \param op_fm_set ...
2138 !> \param mos_new ...
2139 ! **************************************************************************************************
2140  SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
2141 
2142  TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2143  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2144  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2145  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
2146 
2147  INTEGER :: i, icol, lcol, nao, newdim, nmo
2148  LOGICAL :: double_col, double_row
2149  TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
2150  TYPE(cp_fm_type) :: work, work1, work2
2151  TYPE(cp_fm_type), POINTER :: mo_coeff
2152 
2153  DO i = 1, SIZE(op_fm_set, 2) ! spin
2154  CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2155  CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
2156  double_col = .true.
2157  double_row = .false.
2158  CALL cp_fm_struct_double(newstruct, &
2159  mos_new(2*i)%matrix_struct, &
2160  mos_new(2*i)%matrix_struct%context, &
2161  double_col, &
2162  double_row)
2163 
2164  CALL cp_fm_create(work, matrix_struct=newstruct)
2165  CALL cp_fm_create(work1, matrix_struct=newstruct)
2166  CALL cp_fm_create(work2, matrix_struct=newstruct)
2167  CALL cp_fm_get_info(work, ncol_global=newdim)
2168 
2169  CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
2170  DO icol = 1, lcol
2171  work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
2172  work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
2173  END DO
2174 
2175  CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
2176  CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
2177 
2178  DO icol = 1, lcol
2179  work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
2180  work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
2181  END DO
2182 
2183  CALL cp_fm_release(work1)
2184  CALL cp_fm_release(work2)
2185 
2186  CALL cp_fm_struct_double(newstruct1, &
2187  op_fm_set(1, i)%matrix_struct, &
2188  op_fm_set(1, i)%matrix_struct%context, &
2189  double_col, &
2190  double_row)
2191 
2192  CALL cp_fm_create(work1, matrix_struct=newstruct1)
2193 
2194  CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
2195  work, 0.0_dp, work1)
2196 
2197  DO icol = 1, lcol
2198  op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
2199  op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
2200  END DO
2201 
2202  CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
2203  work, 0.0_dp, work1)
2204 
2205  DO icol = 1, lcol
2206  op_fm_set(1, i)%local_data(:, icol) = &
2207  op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
2208  op_fm_set(2, i)%local_data(:, icol) = &
2209  op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
2210  END DO
2211 
2212  CALL cp_fm_release(work)
2213  CALL cp_fm_release(work1)
2214  CALL cp_fm_struct_release(newstruct)
2215  CALL cp_fm_struct_release(newstruct1)
2216 
2217  END DO
2218 
2219  END SUBROUTINE op_orbbas_rtp
2220 
2221 ! **************************************************************************************************
2222 !> \brief ...
2223 !> \param qs_env ...
2224 !> \param magnetic ...
2225 !> \param nmoments ...
2226 !> \param reference ...
2227 !> \param ref_point ...
2228 !> \param unit_number ...
2229 !> \param vel_reprs ...
2230 !> \param com_nl ...
2231 ! **************************************************************************************************
2232  SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
2233 
2234  TYPE(qs_environment_type), POINTER :: qs_env
2235  LOGICAL, INTENT(IN) :: magnetic
2236  INTEGER, INTENT(IN) :: nmoments, reference
2237  REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
2238  INTEGER, INTENT(IN) :: unit_number
2239  LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
2240 
2241  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_locop'
2242 
2243  CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
2244  CHARACTER(LEN=default_string_length) :: description
2245  INTEGER :: akind, handle, i, ia, iatom, idir, &
2246  ikind, ispin, ix, iy, iz, l, nm, nmom, &
2247  order
2248  LOGICAL :: my_com_nl, my_velreprs
2249  REAL(dp) :: charge, dd, strace, trace
2250  REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
2251  nlcom_rv, nlcom_rvr, nlcom_rxrv, &
2252  qupole_der, rmom_vel
2253  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
2254  REAL(dp), DIMENSION(3) :: rcc, ria
2255  TYPE(atomic_kind_type), POINTER :: atomic_kind
2256  TYPE(cell_type), POINTER :: cell
2257  TYPE(cp_result_type), POINTER :: results
2258  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
2259  rho_ao
2260  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
2261  TYPE(dbcsr_type), POINTER :: tmp_ao
2262  TYPE(dft_control_type), POINTER :: dft_control
2263  TYPE(distribution_1d_type), POINTER :: local_particles
2264  TYPE(mp_para_env_type), POINTER :: para_env
2265  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2266  POINTER :: sab_all, sab_orb
2267  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2268  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2269  TYPE(qs_rho_type), POINTER :: rho
2270 
2271  cpassert(ASSOCIATED(qs_env))
2272 
2273  CALL timeset(routinen, handle)
2274 
2275  my_velreprs = .false.
2276  IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
2277  IF (PRESENT(com_nl)) my_com_nl = com_nl
2278  IF (my_velreprs) CALL cite_reference(mattiat2019)
2279 
2280  ! reference point
2281  CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
2282 
2283  ! only allow for moments up to maxl set by basis
2284  nmom = min(nmoments, current_maxl)
2285  ! electronic contribution
2286  NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
2287  CALL get_qs_env(qs_env, &
2288  dft_control=dft_control, &
2289  rho=rho, &
2290  cell=cell, &
2291  results=results, &
2292  particle_set=particle_set, &
2293  qs_kind_set=qs_kind_set, &
2294  para_env=para_env, &
2295  matrix_s=matrix_s, &
2296  sab_all=sab_all, &
2297  sab_orb=sab_orb)
2298 
2299  IF (my_com_nl) THEN
2300  IF ((nmom >= 1) .AND. my_velreprs) THEN
2301  ALLOCATE (nlcom_rv(3))
2302  nlcom_rv(:) = 0._dp
2303  END IF
2304  IF ((nmom >= 2) .AND. my_velreprs) THEN
2305  ALLOCATE (nlcom_rrv(6))
2306  nlcom_rrv(:) = 0._dp
2307  ALLOCATE (nlcom_rvr(6))
2308  nlcom_rvr(:) = 0._dp
2309  ALLOCATE (nlcom_rrv_vrr(6))
2310  nlcom_rrv_vrr(:) = 0._dp
2311  END IF
2312  IF (magnetic) THEN
2313  ALLOCATE (nlcom_rxrv(3))
2314  nlcom_rxrv = 0._dp
2315  END IF
2316  ! Calculate non local correction terms
2317  CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
2318  END IF
2319 
2320  NULLIFY (moments)
2321  nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
2322  CALL dbcsr_allocate_matrix_set(moments, nm)
2323  DO i = 1, nm
2324  ALLOCATE (moments(i)%matrix)
2325  IF (my_velreprs .AND. (nmom >= 2)) THEN
2326  CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
2327  matrix_type=dbcsr_type_symmetric)
2328  CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
2329  ELSE
2330  CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
2331  END IF
2332  CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
2333  END DO
2334 
2335  ! calculate derivatives if quadrupole in vel. reprs. is requested
2336  IF (my_velreprs .AND. (nmom >= 2)) THEN
2337  NULLIFY (moments_der)
2338  CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
2339  DO i = 1, 3 ! x, y, z
2340  DO idir = 1, 3 ! d/dx, d/dy, d/dz
2341  CALL dbcsr_init_p(moments_der(i, idir)%matrix)
2342  CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
2343  matrix_type=dbcsr_type_antisymmetric)
2344  CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
2345  CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
2346  END DO
2347  END DO
2348  CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
2349  ELSE
2350  CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
2351  END IF
2352 
2353  CALL qs_rho_get(rho, rho_ao=rho_ao)
2354 
2355  ALLOCATE (rmom(nm + 1, 3))
2356  ALLOCATE (rlab(nm + 1))
2357  rmom = 0.0_dp
2358  rlab = ""
2359 
2360  IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2361  ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
2362  ! symmetric matrices)
2363  NULLIFY (tmp_ao)
2364  CALL dbcsr_init_p(tmp_ao)
2365  CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2366  CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2367  CALL dbcsr_set(tmp_ao, 0.0_dp)
2368  END IF
2369 
2370  trace = 0.0_dp
2371  DO ispin = 1, dft_control%nspins
2372  CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
2373  rmom(1, 1) = rmom(1, 1) + trace
2374  END DO
2375 
2376  DO i = 1, SIZE(moments)
2377  strace = 0._dp
2378  DO ispin = 1, dft_control%nspins
2379  IF (my_velreprs .AND. nmoments >= 2) THEN
2380  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
2381  0.0_dp, tmp_ao)
2382  CALL dbcsr_trace(tmp_ao, trace)
2383  ELSE
2384  CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
2385  END IF
2386  strace = strace + trace
2387  END DO
2388  rmom(i + 1, 1) = strace
2389  END DO
2390 
2391  CALL dbcsr_deallocate_matrix_set(moments)
2392 
2393  ! nuclear contribution
2394  CALL get_qs_env(qs_env=qs_env, &
2395  local_particles=local_particles)
2396  DO ikind = 1, SIZE(local_particles%n_el)
2397  DO ia = 1, local_particles%n_el(ikind)
2398  iatom = local_particles%list(ikind)%array(ia)
2399  ! fold atomic positions back into unit cell
2400  ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
2401  ria = ria - rcc
2402  atomic_kind => particle_set(iatom)%atomic_kind
2403  CALL get_atomic_kind(atomic_kind, kind_number=akind)
2404  CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
2405  rmom(1, 2) = rmom(1, 2) - charge
2406  DO l = 1, nm
2407  ix = indco(1, l + 1)
2408  iy = indco(2, l + 1)
2409  iz = indco(3, l + 1)
2410  dd = 1._dp
2411  IF (ix > 0) dd = dd*ria(1)**ix
2412  IF (iy > 0) dd = dd*ria(2)**iy
2413  IF (iz > 0) dd = dd*ria(3)**iz
2414  rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2415  CALL set_label(rlab(l + 1), ix, iy, iz)
2416  END DO
2417  END DO
2418  END DO
2419  CALL para_env%sum(rmom(:, 2))
2420  rmom(:, :) = -rmom(:, :)
2421  rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2422 
2423  ! magnetic moments
2424  IF (magnetic) THEN
2425  NULLIFY (magmom)
2426  CALL dbcsr_allocate_matrix_set(magmom, 3)
2427  DO i = 1, SIZE(magmom)
2428  CALL dbcsr_init_p(magmom(i)%matrix)
2429  CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2430  matrix_type=dbcsr_type_antisymmetric)
2431  CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2432  CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2433  END DO
2434 
2435  CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2436 
2437  ALLOCATE (mmom(SIZE(magmom)))
2438  mmom(:) = 0.0_dp
2439  IF (qs_env%run_rtp) THEN
2440  ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2441  ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2442  ! There may be other cases, where the imaginary part of the density is relevant
2443  NULLIFY (rho_ao)
2444  CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2445  END IF
2446  ! if the density is purely real this is an expensive way to calculate zero
2447  DO i = 1, SIZE(magmom)
2448  strace = 0._dp
2449  DO ispin = 1, dft_control%nspins
2450  CALL dbcsr_set(tmp_ao, 0.0_dp)
2451  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2452  0.0_dp, tmp_ao)
2453  CALL dbcsr_trace(tmp_ao, trace)
2454  strace = strace + trace
2455  END DO
2456  mmom(i) = strace
2457  END DO
2458 
2459  CALL dbcsr_deallocate_matrix_set(magmom)
2460  END IF
2461 
2462  ! velocity representations
2463  IF (my_velreprs) THEN
2464  ALLOCATE (rmom_vel(nm))
2465  rmom_vel = 0.0_dp
2466 
2467  DO order = 1, nmom
2468  SELECT CASE (order)
2469 
2470  CASE (1) ! expectation value of momentum
2471  NULLIFY (momentum)
2472  CALL dbcsr_allocate_matrix_set(momentum, 3)
2473  DO i = 1, 3
2474  CALL dbcsr_init_p(momentum(i)%matrix)
2475  CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2476  matrix_type=dbcsr_type_antisymmetric)
2477  CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2478  CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2479  END DO
2480  CALL build_lin_mom_matrix(qs_env, momentum)
2481 
2482  ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2483  IF (qs_env%run_rtp) THEN
2484  NULLIFY (rho_ao)
2485  CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2486  DO idir = 1, SIZE(momentum)
2487  strace = 0._dp
2488  DO ispin = 1, dft_control%nspins
2489  CALL dbcsr_set(tmp_ao, 0.0_dp)
2490  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2491  0.0_dp, tmp_ao)
2492  CALL dbcsr_trace(tmp_ao, trace)
2493  strace = strace + trace
2494  END DO
2495  rmom_vel(idir) = rmom_vel(idir) + strace
2496  END DO
2497  END IF
2498 
2499  CALL dbcsr_deallocate_matrix_set(momentum)
2500 
2501  CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2502  ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2503  qupole_der = 0._dp
2504 
2505  NULLIFY (rho_ao)
2506  CALL qs_rho_get(rho, rho_ao=rho_ao)
2507 
2508  ! Calculate expectation value over real part
2509  trace = 0._dp
2510  DO i = 1, 3
2511  DO idir = 1, 3
2512  strace = 0._dp
2513  DO ispin = 1, dft_control%nspins
2514  CALL dbcsr_set(tmp_ao, 0._dp)
2515  CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2516  CALL dbcsr_trace(tmp_ao, trace)
2517  strace = strace + trace
2518  END DO
2519  qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2520  END DO
2521  END DO
2522 
2523  IF (qs_env%run_rtp) THEN
2524  NULLIFY (rho_ao)
2525  CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2526 
2527  ! Calculate expectation value over imaginary part
2528  trace = 0._dp
2529  DO i = 1, 3
2530  DO idir = 1, 3
2531  strace = 0._dp
2532  DO ispin = 1, dft_control%nspins
2533  CALL dbcsr_set(tmp_ao, 0._dp)
2534  CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2535  CALL dbcsr_trace(tmp_ao, trace)
2536  strace = strace + trace
2537  END DO
2538  qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2539  END DO
2540  END DO
2541  END IF
2542 
2543  ! calculate vel. reprs. of quadrupole moment from derivatives
2544  rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2545  rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2546  rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2547  rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2548  rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2549  rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2550 
2551  DEALLOCATE (qupole_der)
2552  CASE DEFAULT
2553  END SELECT
2554  END DO
2555  END IF
2556 
2557  IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2558  CALL dbcsr_deallocate_matrix(tmp_ao)
2559  END IF
2560  IF (my_velreprs .AND. (nmoments >= 2)) THEN
2561  CALL dbcsr_deallocate_matrix_set(moments_der)
2562  END IF
2563 
2564  description = "[DIPOLE]"
2565  CALL cp_results_erase(results=results, description=description)
2566  CALL put_results(results=results, description=description, &
2567  values=rmom(2:4, 3))
2568 
2569  IF (magnetic .AND. my_velreprs) THEN
2570  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, rmom_vel=rmom_vel)
2571  ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2572  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom)
2573  ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2574  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., rmom_vel=rmom_vel)
2575  ELSE
2576  CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false.)
2577  END IF
2578 
2579  IF (my_com_nl) THEN
2580  IF (magnetic) THEN
2581  mmom(:) = nlcom_rxrv(:)
2582  END IF
2583  IF (my_velreprs .AND. (nmom >= 1)) THEN
2584  DEALLOCATE (rmom_vel)
2585  ALLOCATE (rmom_vel(21))
2586  rmom_vel(1:3) = nlcom_rv
2587  END IF
2588  IF (my_velreprs .AND. (nmom >= 2)) THEN
2589  rmom_vel(4:9) = nlcom_rrv
2590  rmom_vel(10:15) = nlcom_rvr
2591  rmom_vel(16:21) = nlcom_rrv_vrr
2592  END IF
2593  IF (magnetic .AND. .NOT. my_velreprs) THEN
2594  CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2595  ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2596  CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2597  ELSE IF (my_velreprs .AND. magnetic) THEN
2598  CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2599  END IF
2600 
2601  END IF
2602 
2603  IF (my_com_nl) THEN
2604  IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2605  IF (nmom >= 2 .AND. my_velreprs) THEN
2606  DEALLOCATE (nlcom_rrv)
2607  DEALLOCATE (nlcom_rvr)
2608  DEALLOCATE (nlcom_rrv_vrr)
2609  END IF
2610  IF (magnetic) DEALLOCATE (nlcom_rxrv)
2611  END IF
2612 
2613  DEALLOCATE (rmom)
2614  DEALLOCATE (rlab)
2615  IF (magnetic) THEN
2616  DEALLOCATE (mmom)
2617  END IF
2618  IF (my_velreprs) THEN
2619  DEALLOCATE (rmom_vel)
2620  END IF
2621 
2622  CALL timestop(handle)
2623 
2624  END SUBROUTINE qs_moment_locop
2625 
2626 ! **************************************************************************************************
2627 !> \brief ...
2628 !> \param label ...
2629 !> \param ix ...
2630 !> \param iy ...
2631 !> \param iz ...
2632 ! **************************************************************************************************
2633  SUBROUTINE set_label(label, ix, iy, iz)
2634  CHARACTER(LEN=*), INTENT(OUT) :: label
2635  INTEGER, INTENT(IN) :: ix, iy, iz
2636 
2637  INTEGER :: i
2638 
2639  label = ""
2640  DO i = 1, ix
2641  WRITE (label(i:), "(A1)") "X"
2642  END DO
2643  DO i = ix + 1, ix + iy
2644  WRITE (label(i:), "(A1)") "Y"
2645  END DO
2646  DO i = ix + iy + 1, ix + iy + iz
2647  WRITE (label(i:), "(A1)") "Z"
2648  END DO
2649 
2650  END SUBROUTINE set_label
2651 
2652 ! **************************************************************************************************
2653 !> \brief ...
2654 !> \param unit_number ...
2655 !> \param nmom ...
2656 !> \param rmom ...
2657 !> \param rlab ...
2658 !> \param rcc ...
2659 !> \param cell ...
2660 !> \param periodic ...
2661 !> \param mmom ...
2662 !> \param rmom_vel ...
2663 ! **************************************************************************************************
2664  SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2665  INTEGER, INTENT(IN) :: unit_number, nmom
2666  REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2667  CHARACTER(LEN=8), DIMENSION(:) :: rlab
2668  REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2669  TYPE(cell_type), POINTER :: cell
2670  LOGICAL :: periodic
2671  REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2672 
2673  INTEGER :: i, i0, i1, j, l
2674  REAL(dp) :: dd
2675 
2676  IF (unit_number > 0) THEN
2677  DO l = 0, nmom
2678  SELECT CASE (l)
2679  CASE (0)
2680  WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2681  WRITE (unit_number, "(T3,A)") "Charges"
2682  WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2683  "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2684  CASE (1)
2685  IF (periodic) THEN
2686  WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2687  WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2688  WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2689  WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2690  WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2691  ELSE
2692  WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2693  END IF
2694  dd = sqrt(sum(rmom(2:4, 3)**2))*debye
2695  WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2696  WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2697  (trim(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2698  CASE (2)
2699  WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2700  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2701  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2702  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2703  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2704  CASE (3)
2705  WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2706  WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2707  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2708  WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2709  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2710  WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2711  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2712  CASE (4)
2713  WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2714  WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2715  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2716  WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2717  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2718  WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2719  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2720  WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2721  (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2722  CASE DEFAULT
2723  WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2724  " L=", l
2725  i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2726  i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2727  dd = debye/(bohr)**(l - 1)
2728  DO i = i0, i1, 3
2729  WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2730  (trim(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, min(i1, i + 2))
2731  END DO
2732  END SELECT
2733  END DO
2734  IF (PRESENT(mmom)) THEN
2735  IF (nmom >= 1) THEN
2736  dd = sqrt(sum(mmom(1:3)**2))
2737  WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2738  WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2739  (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2740  END IF
2741  END IF
2742  IF (PRESENT(rmom_vel)) THEN
2743  DO l = 1, nmom
2744  SELECT CASE (l)
2745  CASE (1)
2746  dd = sqrt(sum(rmom_vel(1:3)**2))
2747  WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2748  WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2749  (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2750  CASE (2)
2751  WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2752  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2753  (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2754  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2755  (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2756  CASE DEFAULT
2757  END SELECT
2758  END DO
2759  END IF
2760  END IF
2761 
2762  END SUBROUTINE print_moments
2763 
2764 ! **************************************************************************************************
2765 !> \brief ...
2766 !> \param unit_number ...
2767 !> \param nmom ...
2768 !> \param rlab ...
2769 !> \param mmom ...
2770 !> \param rmom_vel ...
2771 ! **************************************************************************************************
2772  SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2773  INTEGER, INTENT(IN) :: unit_number, nmom
2774  CHARACTER(LEN=8), DIMENSION(:) :: rlab
2775  REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2776 
2777  INTEGER :: i, l
2778  REAL(dp) :: dd
2779 
2780  IF (unit_number > 0) THEN
2781  IF (PRESENT(mmom)) THEN
2782  IF (nmom >= 1) THEN
2783  dd = sqrt(sum(mmom(1:3)**2))
2784  WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2785  WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2786  (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2787  END IF
2788  END IF
2789  IF (PRESENT(rmom_vel)) THEN
2790  DO l = 1, nmom
2791  SELECT CASE (l)
2792  CASE (1)
2793  dd = sqrt(sum(rmom_vel(1:3)**2))
2794  WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2795  WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2796  (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2797  CASE (2)
2798  WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2799  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2800  (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2801  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2802  (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2803  WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2804  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2805  (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2806  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2807  (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2808  WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2809  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2810  (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2811  WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2812  (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2813  CASE DEFAULT
2814  END SELECT
2815  END DO
2816  END IF
2817  END IF
2818 
2819  END SUBROUTINE print_moments_nl
2820 
2821 ! **************************************************************************************************
2822 !> \brief Calculate the expectation value of operators related to non-local potential:
2823 !> [r, Vnl], noted rv
2824 !> r x [r,Vnl], noted rxrv
2825 !> [rr,Vnl], noted rrv
2826 !> r x Vnl x r, noted rvr
2827 !> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2828 !> Note that the 3 first operator are commutator while the 2 last
2829 !> are not. For reading clarity the same notation is used for all 5
2830 !> operators.
2831 !> \param qs_env ...
2832 !> \param nlcom_rv ...
2833 !> \param nlcom_rxrv ...
2834 !> \param nlcom_rrv ...
2835 !> \param nlcom_rvr ...
2836 !> \param nlcom_rrv_vrr ...
2837 !> \param ref_point ...
2838 ! **************************************************************************************************
2839  SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2840  nlcom_rrv_vrr, ref_point)
2841 
2842  TYPE(qs_environment_type), POINTER :: qs_env
2843  REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2844  nlcom_rvr, nlcom_rrv_vrr
2845  REAL(dp), DIMENSION(3) :: ref_point
2846 
2847  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_commutator_nl_terms'
2848 
2849  INTEGER :: handle, ind, ispin
2850  LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2851  calc_rvr, calc_rxrv
2852  REAL(dp) :: eps_ppnl, strace, trace
2853  TYPE(cell_type), POINTER :: cell
2854  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2855  matrix_rvr, matrix_rxrv, matrix_s, &
2856  rho_ao
2857  TYPE(dbcsr_type), POINTER :: tmp_ao
2858  TYPE(dft_control_type), POINTER :: dft_control
2859  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2860  POINTER :: sab_all, sab_orb, sap_ppnl
2861  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2862  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2863  TYPE(qs_rho_type), POINTER :: rho
2864 
2865  CALL timeset(routinen, handle)
2866 
2867  calc_rv = .false.
2868  calc_rxrv = .false.
2869  calc_rrv = .false.
2870  calc_rvr = .false.
2871  calc_rrv_vrr = .false.
2872 
2873  ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2874  ! The real part of the density matrix rho_ao is symmetric so that
2875  ! the expectation value of real density matrix is zero. Hence, if
2876  ! the density matrix is real, no need to compute these quantities.
2877  ! This is not the case for rvr and rrv_vrr which are symmetric.
2878 
2879  IF (ALLOCATED(nlcom_rv)) THEN
2880  nlcom_rv(:) = 0._dp
2881  IF (qs_env%run_rtp) calc_rv = .true.
2882  END IF
2883  IF (ALLOCATED(nlcom_rxrv)) THEN
2884  nlcom_rxrv(:) = 0._dp
2885  IF (qs_env%run_rtp) calc_rxrv = .true.
2886  END IF
2887  IF (ALLOCATED(nlcom_rrv)) THEN
2888  nlcom_rrv(:) = 0._dp
2889  IF (qs_env%run_rtp) calc_rrv = .true.
2890  END IF
2891  IF (ALLOCATED(nlcom_rvr)) THEN
2892  nlcom_rvr(:) = 0._dp
2893  calc_rvr = .true.
2894  END IF
2895  IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2896  nlcom_rrv_vrr(:) = 0._dp
2897  calc_rrv_vrr = .true.
2898  END IF
2899 
2900  IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv &
2901  .OR. calc_rvr .OR. calc_rrv_vrr)) RETURN
2902 
2903  NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2904  CALL get_qs_env(qs_env, &
2905  cell=cell, &
2906  dft_control=dft_control, &
2907  matrix_s=matrix_s, &
2908  particle_set=particle_set, &
2909  qs_kind_set=qs_kind_set, &
2910  rho=rho, &
2911  sab_orb=sab_orb, &
2912  sab_all=sab_all, &
2913  sap_ppnl=sap_ppnl)
2914 
2915  eps_ppnl = dft_control%qs_control%eps_ppnl
2916 
2917  ! Allocate storage
2918  NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2919  IF (calc_rv) THEN
2920  CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2921  DO ind = 1, 3
2922  CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2923  CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2924  matrix_type=dbcsr_type_antisymmetric)
2925  CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2926  CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2927  END DO
2928  END IF
2929 
2930  IF (calc_rxrv) THEN
2931  CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2932  DO ind = 1, 3
2933  CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2934  CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2935  matrix_type=dbcsr_type_antisymmetric)
2936  CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2937  CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2938  END DO
2939  END IF
2940 
2941  IF (calc_rrv) THEN
2942  CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2943  DO ind = 1, 6
2944  CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2945  CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2946  matrix_type=dbcsr_type_antisymmetric)
2947  CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2948  CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2949  END DO
2950  END IF
2951 
2952  IF (calc_rvr) THEN
2953  CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2954  DO ind = 1, 6
2955  CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2956  CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2957  matrix_type=dbcsr_type_symmetric)
2958  CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2959  CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2960  END DO
2961  END IF
2962  IF (calc_rrv_vrr) THEN
2963  CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2964  DO ind = 1, 6
2965  CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2966  CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2967  matrix_type=dbcsr_type_symmetric)
2968  CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2969  CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2970  END DO
2971  END IF
2972 
2973  ! calculate evaluation of operators in AO basis set
2974  CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2975  matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2976  matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2977 
2978  ! Calculate expectation values
2979  ! Real part
2980  NULLIFY (tmp_ao)
2981  CALL dbcsr_init_p(tmp_ao)
2982  CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2983  CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2984  CALL dbcsr_set(tmp_ao, 0.0_dp)
2985 
2986  IF (calc_rvr .OR. calc_rrv_vrr) THEN
2987  NULLIFY (rho_ao)
2988  CALL qs_rho_get(rho, rho_ao=rho_ao)
2989 
2990  IF (calc_rvr) THEN
2991  trace = 0._dp
2992  DO ind = 1, SIZE(matrix_rvr)
2993  strace = 0._dp
2994  DO ispin = 1, dft_control%nspins
2995  CALL dbcsr_set(tmp_ao, 0.0_dp)
2996  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2997  0.0_dp, tmp_ao)
2998  CALL dbcsr_trace(tmp_ao, trace)
2999  strace = strace + trace
3000  END DO
3001  nlcom_rvr(ind) = nlcom_rvr(ind) + strace
3002  END DO
3003  END IF
3004 
3005  IF (calc_rrv_vrr) THEN
3006  trace = 0._dp
3007  DO ind = 1, SIZE(matrix_rrv_vrr)
3008  strace = 0._dp
3009  DO ispin = 1, dft_control%nspins
3010  CALL dbcsr_set(tmp_ao, 0.0_dp)
3011  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
3012  0.0_dp, tmp_ao)
3013  CALL dbcsr_trace(tmp_ao, trace)
3014  strace = strace + trace
3015  END DO
3016  nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
3017  END DO
3018  END IF
3019  END IF
3020 
3021  ! imagninary part of the density matrix
3022  NULLIFY (rho_ao)
3023  CALL qs_rho_get(rho, rho_ao_im=rho_ao)
3024 
3025  IF (calc_rv) THEN
3026  trace = 0._dp
3027  DO ind = 1, SIZE(matrix_rv)
3028  strace = 0._dp
3029  DO ispin = 1, dft_control%nspins
3030  CALL dbcsr_set(tmp_ao, 0.0_dp)
3031  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
3032  0.0_dp, tmp_ao)
3033  CALL dbcsr_trace(tmp_ao, trace)
3034  strace = strace + trace
3035  END DO
3036  nlcom_rv(ind) = nlcom_rv(ind) + strace
3037  END DO
3038  END IF
3039 
3040  IF (calc_rrv) THEN
3041  trace = 0._dp
3042  DO ind = 1, SIZE(matrix_rrv)
3043  strace = 0._dp
3044  DO ispin = 1, dft_control%nspins
3045  CALL dbcsr_set(tmp_ao, 0.0_dp)
3046  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
3047  0.0_dp, tmp_ao)
3048  CALL dbcsr_trace(tmp_ao, trace)
3049  strace = strace + trace
3050  END DO
3051  nlcom_rrv(ind) = nlcom_rrv(ind) + strace
3052  END DO
3053  END IF
3054 
3055  IF (calc_rxrv) THEN
3056  trace = 0._dp
3057  DO ind = 1, SIZE(matrix_rxrv)
3058  strace = 0._dp
3059  DO ispin = 1, dft_control%nspins
3060  CALL dbcsr_set(tmp_ao, 0.0_dp)
3061  CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
3062  0.0_dp, tmp_ao)
3063  CALL dbcsr_trace(tmp_ao, trace)
3064  strace = strace + trace
3065  END DO
3066  nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
3067  END DO
3068  END IF
3069  CALL dbcsr_deallocate_matrix(tmp_ao)
3070  IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
3071  IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
3072  IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
3073  IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
3074  IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
3075 
3076  CALL timestop(handle)
3077  END SUBROUTINE calculate_commutator_nl_terms
3078 
3079 ! *****************************************************************************
3080 !> \brief ...
3081 !> \param qs_env ...
3082 !> \param difdip ...
3083 !> \param deltaR ...
3084 !> \param order ...
3085 !> \param rcc ...
3086 !> \note calculate matrix elements <a|r_beta |db/dR_alpha > + <da/dR_alpha | r_beta | b >
3087 !> be aware: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > only valid
3088 !> if alpha .neq.beta
3089 !> if alpha=beta: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > - < a | b >
3090 !> modified from qs_efield_mo_derivatives
3091 !> SL July 2015
3092 ! **************************************************************************************************
3093  SUBROUTINE dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
3094  TYPE(qs_environment_type), POINTER :: qs_env
3095  TYPE(dbcsr_p_type), DIMENSION(:, :), &
3096  INTENT(INOUT), POINTER :: difdip
3097  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
3098  POINTER :: deltar
3099  INTEGER, INTENT(IN) :: order
3100  REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3101 
3102  CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_deriv_ao'
3103 
3104  INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
3105  last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
3106  sgfb
3107  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3108  npgfb, nsgfa, nsgfb
3109  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3110  LOGICAL :: found
3111  REAL(dp) :: dab
3112  REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3113  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3114  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
3115  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
3116  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
3117  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
3118  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: difmab2
3119  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
3120  TYPE(cell_type), POINTER :: cell
3121  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3122  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3123  TYPE(neighbor_list_iterator_p_type), &
3124  DIMENSION(:), POINTER :: nl_iterator
3125  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3126  POINTER :: sab_all, sab_orb
3127  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3128  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3129  TYPE(qs_kind_type), POINTER :: qs_kind
3130 
3131  CALL timeset(routinen, handle)
3132 
3133  NULLIFY (cell, particle_set, qs_kind_set, sab_orb, sab_all)
3134  CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
3135  qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_all=sab_all)
3136  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
3137  maxco=ldab, maxsgf=maxsgf)
3138 
3139  nkind = SIZE(qs_kind_set)
3140  natom = SIZE(particle_set)
3141 
3142  m_dim = ncoset(order) - 1
3143 
3144  IF (PRESENT(rcc)) THEN
3145  rc = rcc
3146  ELSE
3147  rc = 0._dp
3148  END IF
3149 
3150  ALLOCATE (basis_set_list(nkind))
3151 
3152  ALLOCATE (mab(ldab, ldab, m_dim))
3153  ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
3154  ALLOCATE (work(ldab, maxsgf))
3155  ALLOCATE (mint(3, 3))
3156  ALLOCATE (mint2(3, 3))
3157 
3158  mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
3159  difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
3160  work(1:ldab, 1:maxsgf) = 0.0_dp
3161 
3162  DO i = 1, 3
3163  DO j = 1, 3
3164  NULLIFY (mint(i, j)%block)
3165  NULLIFY (mint2(i, j)%block)
3166  END DO
3167  END DO
3168 
3169  ! Set the basis_set_list(nkind) to point to the corresponding basis sets
3170  DO ikind = 1, nkind
3171  qs_kind => qs_kind_set(ikind)
3172  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
3173  IF (ASSOCIATED(basis_set_a)) THEN
3174  basis_set_list(ikind)%gto_basis_set => basis_set_a
3175  ELSE
3176  NULLIFY (basis_set_list(ikind)%gto_basis_set)
3177  END IF
3178  END DO
3179 
3180  CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3181  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3182  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
3183  iatom=iatom, jatom=jatom, r=rab)
3184 
3185  basis_set_a => basis_set_list(ikind)%gto_basis_set
3186  basis_set_b => basis_set_list(jkind)%gto_basis_set
3187  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3188  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3189 
3190  ! basis ikind
3191  first_sgfa => basis_set_a%first_sgf
3192  la_max => basis_set_a%lmax
3193  la_min => basis_set_a%lmin
3194  npgfa => basis_set_a%npgf
3195  nseta = basis_set_a%nset
3196  nsgfa => basis_set_a%nsgf_set
3197  rpgfa => basis_set_a%pgf_radius
3198  set_radius_a => basis_set_a%set_radius
3199  sphi_a => basis_set_a%sphi
3200  zeta => basis_set_a%zet
3201  ! basis jkind
3202  first_sgfb => basis_set_b%first_sgf
3203  lb_max => basis_set_b%lmax
3204  lb_min => basis_set_b%lmin
3205  npgfb => basis_set_b%npgf
3206  nsetb = basis_set_b%nset
3207  nsgfb => basis_set_b%nsgf_set
3208  rpgfb => basis_set_b%pgf_radius
3209  set_radius_b => basis_set_b%set_radius
3210  sphi_b => basis_set_b%sphi
3211  zetb => basis_set_b%zet
3212 
3213  IF (inode == 1) last_jatom = 0
3214 
3215  ! this guarentees minimum image convention
3216  ! anything else would not make sense
3217  IF (jatom == last_jatom) THEN
3218  cycle
3219  END IF
3220 
3221  last_jatom = jatom
3222 
3223  irow = iatom
3224  icol = jatom
3225 
3226  DO i = 1, 3
3227  DO j = 1, 3
3228  NULLIFY (mint(i, j)%block)
3229  CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
3230  row=irow, col=icol, block=mint(i, j)%block, &
3231  found=found)
3232  cpassert(found)
3233  END DO
3234  END DO
3235 
3236  ra(:) = particle_set(iatom)%r(:)
3237  rb(:) = particle_set(jatom)%r(:)
3238  rab(:) = pbc(rb, ra, cell)
3239  rac(:) = pbc(ra - rc, cell)
3240  rbc(:) = pbc(rb - rc, cell)
3241  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
3242 
3243  DO iset = 1, nseta
3244  ncoa = npgfa(iset)*ncoset(la_max(iset))
3245  sgfa = first_sgfa(1, iset)
3246  DO jset = 1, nsetb
3247  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
3248  ncob = npgfb(jset)*ncoset(lb_max(jset))
3249  sgfb = first_sgfb(1, jset)
3250  ldab = max(ncoa, ncob)
3251  lda = ncoset(la_max(iset))*npgfa(iset)
3252  ldb = ncoset(lb_max(jset))*npgfb(jset)
3253  ALLOCATE (difmab(lda, ldb, m_dim, 3))
3254 
3255  ! Calculate integral (da|r|b)
3256  CALL diff_momop2(la_max(iset), npgfa(iset), zeta(:, iset), &
3257  rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
3258  zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
3259  difmab, deltar=deltar, iatom=iatom, jatom=jatom)
3260 
3261 ! *** Contraction step ***
3262 
3263  DO idir = 1, 3 ! derivative of AO function
3264  DO j = 1, 3 ! position operator r_j
3265  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3266  1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
3267  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
3268  0.0_dp, work(1, 1), SIZE(work, 1))
3269 
3270  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3271  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
3272  work(1, 1), SIZE(work, 1), &
3273  1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
3274  SIZE(mint(j, idir)%block, 1))
3275  END DO !j
3276  END DO !idir
3277  DEALLOCATE (difmab)
3278  END DO !jset
3279  END DO !iset
3280  END do!iterator
3281 
3282  CALL neighbor_list_iterator_release(nl_iterator)
3283 
3284  DO i = 1, 3
3285  DO j = 1, 3
3286  NULLIFY (mint(i, j)%block)
3287  END DO
3288  END DO
3289 
3290  DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
3291 
3292  CALL timestop(handle)
3293  END SUBROUTINE dipole_deriv_ao
3294 
3295 END MODULE qs_moments
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
ncoset
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition: grid_common.h:73
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
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition: ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
...
Definition: ai_angmom.F:52
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition: ai_moments.F:17
ai_moments::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::diff_momop_velocity
subroutine, public diff_momop_velocity(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, lambda, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the primitive on the r...
Definition: ai_moments.F:79
ai_moments::contract_cossin
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition: ai_moments.F:261
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition: ai_moments.F:339
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_moments::diff_momop2
subroutine, public diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition: ai_moments.F:164
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::mattiat2019
integer, save, public mattiat2019
Definition: bibliography.F:43
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
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition: commutator_rpnl.F:439
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition: cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_det
subroutine, public cp_cfm_det(matrix_a, det_a)
Computes the determinant (with a correct sign even in parallel environment!) of a complex square matr...
Definition: cp_cfm_basic_linalg.F:76
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_double
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
Definition: cp_fm_struct.F:536
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_methods.F:17
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition: cp_result_methods.F:362
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition: cp_result_types.F:18
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
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
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition: moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition: moments_utils.F:61
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::current_maxl
integer, save, public current_maxl
Definition: orbital_pointers.F:40
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
physcon::debye
real(kind=dp), parameter, public debye
Definition: physcon.F:201
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_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
qs_moments::build_local_magmom_matrix
subroutine, public build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:1123
qs_moments::dipole_velocity_deriv
subroutine, public dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the basis function on ...
Definition: qs_moments.F:126
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:558
qs_moments::qs_moment_locop
subroutine, public qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
...
Definition: qs_moments.F:2233
qs_moments::dipole_deriv_ao
subroutine, public dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
...
Definition: qs_moments.F:3094
qs_moments::build_berry_moment_matrix
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Definition: qs_moments.F:1336
qs_moments::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition: qs_moments.F:790
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition: qs_moments.F:1503
qs_moments::build_dsdv_moments
subroutine, public build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
Builds the moments for the derivative of the overlap with respect to nuclear velocities.
Definition: qs_moments.F:329
qs_moments::qs_moment_berry_phase
subroutine, public qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
...
Definition: qs_moments.F:1712
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_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_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition: rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_H_old, exp_H_new, H_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, S_inv, S_half, S_minus_half, B_mat, C_mat, propagator_matrix, mixing, mixing_factor, S_der, dt, nsteps, SinvH, SinvH_imag, SinvB, admm_mos)
...
Definition: rt_propagation_types.F:408