LCOV - code coverage report
Current view: top level - src - qs_moments.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:591cf04) Lines: 93.0 % 1623 1510
Test Date: 2026-09-21 02:17:57 Functions: 100.0 % 18 18

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \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              :                          moment
      20              :    USE atomic_kind_types, ONLY: atomic_kind_type, &
      21              :                                 get_atomic_kind
      22              :    USE basis_set_types, ONLY: gto_basis_set_p_type, &
      23              :                               gto_basis_set_type
      24              :    USE bibliography, ONLY: Mattiat2019, &
      25              :                            cite_reference
      26              :    USE block_p_types, ONLY: block_p_type
      27              :    USE cell_types, ONLY: cell_type, &
      28              :                          pbc, &
      29              :                          get_cell
      30              :    USE commutator_rpnl, ONLY: build_com_mom_nl
      31              :    USE cp_blacs_env, ONLY: cp_blacs_env_type
      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_api, ONLY: &
      39              :       dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
      40              :       dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
      41              :       dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
      42              :    USE cp_dbcsr_contrib, ONLY: dbcsr_dot, &
      43              :                                dbcsr_trace
      44              :    USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
      45              :    USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply, &
      46              :                                   dbcsr_allocate_matrix_set, &
      47              :                                   dbcsr_deallocate_matrix_set
      48              :    USE cp_fm_struct, ONLY: cp_fm_struct_create, &
      49              :                            cp_fm_struct_double, &
      50              :                            cp_fm_struct_release, &
      51              :                            cp_fm_struct_type
      52              :    USE cp_fm_types, ONLY: cp_fm_copy_general, &
      53              :                           cp_fm_create, &
      54              :                           cp_fm_get_element, &
      55              :                           cp_fm_get_info, &
      56              :                           cp_fm_release, &
      57              :                           cp_fm_set_all, &
      58              :                           cp_fm_type
      59              :    USE cp_result_methods, ONLY: cp_results_erase, &
      60              :                                 put_results
      61              :    USE cp_result_types, ONLY: cp_result_type
      62              :    USE distribution_1d_types, ONLY: distribution_1d_type
      63              :    USE kinds, ONLY: default_string_length, &
      64              :                     dp, &
      65              :                     max_line_length
      66              :    USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices, &
      67              :                                       rs_to_kp
      68              :    USE kpoint_methods, ONLY: kpoint_init_cell_index, &
      69              :                              kpoint_initialize, &
      70              :                              rskp_transform
      71              :    USE kpoint_types, ONLY: get_kpoint_info, &
      72              :                            kpoint_env_p_type, &
      73              :                            kpoint_env_type, &
      74              :                            kpoint_type, &
      75              :                            kpoint_create, &
      76              :                            kpoint_release, &
      77              :                            read_kpoint_section
      78              :    USE mathconstants, ONLY: pi, &
      79              :                             twopi, &
      80              :                             gaussi, &
      81              :                             z_zero
      82              :    USE message_passing, ONLY: mp_para_env_type
      83              :    USE moments_utils, ONLY: get_reference_point
      84              :    USE orbital_pointers, ONLY: current_maxl, &
      85              :                                indco, &
      86              :                                ncoset
      87              :    USE parallel_gemm_api, ONLY: parallel_gemm
      88              :    USE particle_methods, ONLY: get_particle_set
      89              :    USE particle_types, ONLY: particle_type
      90              :    USE physcon, ONLY: bohr, &
      91              :                       debye, &
      92              :                       angstrom
      93              :    USE qs_environment_types, ONLY: get_qs_env, &
      94              :                                    qs_environment_type
      95              :    USE qs_kind_types, ONLY: get_qs_kind, &
      96              :                             get_qs_kind_set, &
      97              :                             qs_kind_type
      98              :    USE qs_ks_types, ONLY: get_ks_env, &
      99              :                           qs_ks_env_type
     100              :    USE qs_mo_types, ONLY: get_mo_set, &
     101              :                           mo_set_type
     102              :    USE qs_neighbor_list_types, ONLY: get_iterator_info, &
     103              :                                      neighbor_list_iterate, &
     104              :                                      neighbor_list_iterator_create, &
     105              :                                      neighbor_list_iterator_p_type, &
     106              :                                      neighbor_list_iterator_release, &
     107              :                                      neighbor_list_set_p_type
     108              :    USE qs_operators_ao, ONLY: build_exp_ikr_matrix, &
     109              :                               build_lin_mom_matrix
     110              :    USE qs_overlap, ONLY: build_overlap_matrix
     111              :    USE qs_rho_types, ONLY: qs_rho_get, &
     112              :                            qs_rho_type
     113              :    USE rt_propagation_types, ONLY: get_rtp, &
     114              :                                    rt_prop_type
     115              :    USE cp_parser_methods, ONLY: read_float_object
     116              :    USE input_section_types, ONLY: section_vals_get, &
     117              :                                   section_vals_get_subs_vals, &
     118              :                                   section_vals_type, &
     119              :                                   section_vals_val_get
     120              :    USE mathlib, ONLY: geeig_right, &
     121              :                       gemm_square
     122              :    USE string_utilities, ONLY: uppercase
     123              : 
     124              : #include "./base/base_uses.f90"
     125              : 
     126              :    IMPLICIT NONE
     127              : 
     128              :    PRIVATE
     129              : 
     130              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
     131              : 
     132              :    ! Public subroutines
     133              :    PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
     134              :    PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
     135              :    PUBLIC :: qs_moment_berry_phase, qs_moment_locop
     136              :    PUBLIC :: qs_moment_kpoints, build_local_moment_matrix_rs_img
     137              :    PUBLIC :: qs_moment_kpoints_deep
     138              :    PUBLIC :: qs_moment_kpoints_scf_mos
     139              :    PUBLIC :: build_local_moments_der_matrix
     140              :    PUBLIC :: calculate_commutator_nl_terms, op_orbbas, op_orbbas_rtp, &
     141              :              print_moments, print_moments_nl, set_label
     142              : 
     143              : CONTAINS
     144              : 
     145              : ! **************************************************************************************************
     146              : !> \brief ...
     147              : !> \param qs_env ...
     148              : !> \param moments ...
     149              : !> \param nmoments ...
     150              : !> \param ref_point ...
     151              : !> \param ref_points ...
     152              : !> \param basis_type ...
     153              : !> \param all_images accumulate all neighbor images belonging to one atom pair
     154              : !> \param minimum_image apply the half-cell image filter
     155              : !> \param neighbor_image use the neighbor-list image vector for one atom-pair representative
     156              : !> \param first_component first Cartesian moment component to write
     157              : ! **************************************************************************************************
     158         3428 :    SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, &
     159              :                                         all_images, minimum_image, neighbor_image, first_component)
     160              : 
     161              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     162              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: moments
     163              :       INTEGER, INTENT(IN)                                :: nmoments
     164              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: ref_point
     165              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     166              :          OPTIONAL                                        :: ref_points
     167              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
     168              :       LOGICAL, OPTIONAL                                   :: all_images, minimum_image, neighbor_image
     169              :       INTEGER, OPTIONAL                                   :: first_component
     170              : 
     171              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moment_matrix'
     172              : 
     173              :       INTEGER :: component_start, handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, &
     174              :                  last_jatom, maxco, maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
     175              :       LOGICAL                                            :: found, my_all_images, my_minimum_image, &
     176              :                                                             my_neighbor_image, new_atom_pair
     177              :       REAL(KIND=dp)                                      :: dab
     178              :       REAL(KIND=dp), DIMENSION(3)                        :: half_cell
     179         3428 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     180         3428 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: mab
     181              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, rb, rbc, rc
     182         3428 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: mint
     183              :       TYPE(cell_type), POINTER                           :: cell
     184         3428 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     185              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     186              :       TYPE(neighbor_list_iterator_p_type), &
     187         3428 :          DIMENSION(:), POINTER                           :: nl_iterator
     188              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     189         3428 :          POINTER                                         :: sab_orb
     190         3428 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     191         3428 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     192              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     193              : 
     194         3428 :       IF (nmoments < 1) RETURN
     195              : 
     196         3428 :       CALL timeset(routineN, handle)
     197              : 
     198         3428 :       my_all_images = .FALSE.
     199         3428 :       IF (PRESENT(all_images)) my_all_images = all_images
     200         3428 :       my_minimum_image = .FALSE.
     201         3428 :       IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
     202         3428 :       my_neighbor_image = .FALSE.
     203         3428 :       IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
     204            6 :       CPASSERT(.NOT. (my_all_images .AND. my_neighbor_image))
     205              : 
     206         3428 :       nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
     207         3428 :       CPASSERT(SIZE(moments) >= nm)
     208         3428 :       component_start = 1
     209         3428 :       IF (PRESENT(first_component)) component_start = first_component
     210         3428 :       CPASSERT(component_start >= 1 .AND. component_start <= nm)
     211              : 
     212         3428 :       NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
     213              :       CALL get_qs_env(qs_env=qs_env, &
     214              :                       qs_kind_set=qs_kind_set, &
     215              :                       particle_set=particle_set, cell=cell, &
     216         3428 :                       sab_orb=sab_orb)
     217              : 
     218         3428 :       IF (my_minimum_image) THEN
     219          150 :          half_cell = 0.5_dp*NORM2(cell%hmat, DIM=1)
     220              :       END IF
     221              : 
     222         3428 :       nkind = SIZE(qs_kind_set)
     223              :       ! Allocate work storage
     224              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     225              :                            maxco=maxco, maxsgf=maxsgf, &
     226         6750 :                            basis_type=basis_type)
     227              : 
     228        17140 :       ALLOCATE (mab(maxco, maxco, nm))
     229         3428 :       mab(:, :, :) = 0.0_dp
     230              : 
     231        13712 :       ALLOCATE (work(maxco, maxsgf))
     232         3428 :       work(:, :) = 0.0_dp
     233              : 
     234        20738 :       ALLOCATE (mint(nm))
     235        13882 :       DO i = component_start, nm
     236        13882 :          NULLIFY (mint(i)%block)
     237              :       END DO
     238              : 
     239        17166 :       ALLOCATE (basis_set_list(nkind))
     240        10310 :       DO ikind = 1, nkind
     241         6882 :          qs_kind => qs_kind_set(ikind)
     242         6882 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
     243        10310 :          IF (ASSOCIATED(basis_set_a)) THEN
     244         6882 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     245              :          ELSE
     246            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     247              :          END IF
     248              :       END DO
     249         3428 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     250        39522 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     251              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
     252        36094 :                                 iatom=iatom, jatom=jatom, r=rab)
     253        36094 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     254        36094 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     255        36094 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     256        36094 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     257              :          ASSOCIATE ( &
     258              :             ! basis ikind
     259              :             first_sgfa => basis_set_a%first_sgf, &
     260              :             la_max => basis_set_a%lmax, &
     261              :             la_min => basis_set_a%lmin, &
     262              :             npgfa => basis_set_a%npgf, &
     263              :             nsgfa => basis_set_a%nsgf_set, &
     264              :             rpgfa => basis_set_a%pgf_radius, &
     265              :             set_radius_a => basis_set_a%set_radius, &
     266              :             sphi_a => basis_set_a%sphi, &
     267              :             zeta => basis_set_a%zet, &
     268              :             ! basis jkind, &
     269              :             first_sgfb => basis_set_b%first_sgf, &
     270              :             lb_max => basis_set_b%lmax, &
     271              :             lb_min => basis_set_b%lmin, &
     272              :             npgfb => basis_set_b%npgf, &
     273              :             nsgfb => basis_set_b%nsgf_set, &
     274              :             rpgfb => basis_set_b%pgf_radius, &
     275              :             set_radius_b => basis_set_b%set_radius, &
     276              :             sphi_b => basis_set_b%sphi, &
     277              :             zetb => basis_set_b%zet)
     278              : 
     279        36094 :             nseta = basis_set_a%nset
     280        36094 :             nsetb = basis_set_b%nset
     281              : 
     282         9506 :             IF (inode == 1) last_jatom = 0
     283              : 
     284        36094 :             IF (my_minimum_image) THEN
     285          482 :                IF (ANY(ABS(rab(:)) > half_cell(:))) CYCLE
     286              :             END IF
     287              : 
     288        36092 :             new_atom_pair = jatom /= last_jatom
     289        36092 :             IF (.NOT. my_all_images .AND. .NOT. new_atom_pair) THEN
     290              :                CYCLE
     291              :             END IF
     292              : 
     293        13612 :             IF (new_atom_pair) THEN
     294        13363 :                last_jatom = jatom
     295              : 
     296        13363 :                IF (iatom <= jatom) THEN
     297         8323 :                   irow = iatom
     298         8323 :                   icol = jatom
     299              :                ELSE
     300         5040 :                   irow = jatom
     301         5040 :                   icol = iatom
     302              :                END IF
     303              : 
     304        54322 :                DO i = component_start, nm
     305        40959 :                   NULLIFY (mint(i)%block)
     306              :                   CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
     307        40959 :                                          row=irow, col=icol, BLOCK=mint(i)%block, found=found)
     308        40959 :                   CPASSERT(found .AND. ASSOCIATED(mint(i)%block))
     309      2009108 :                   mint(i)%block = 0._dp
     310              :                END DO
     311              :             END IF
     312              : 
     313              :             ! Reference point and image geometry.
     314        13612 :             IF (PRESENT(ref_points)) THEN
     315            0 :                rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
     316        13612 :             ELSE IF (PRESENT(ref_point)) THEN
     317        51368 :                rc(:) = ref_point(:)
     318              :             ELSE
     319          770 :                rc(:) = 0._dp
     320              :             END IF
     321              : 
     322        13612 :             IF (my_all_images) THEN
     323              :                ! Keep the image represented by the current neighbor-list node.
     324          400 :                ra(:) = pbc(particle_set(iatom)%r(:), cell)
     325         1600 :                rb(:) = ra(:) + rab(:)
     326          400 :                rac(:) = pbc(rc, ra, cell)
     327          400 :                rbc(:) = pbc(rc, rb, cell)
     328         1600 :                dab = NORM2(rab)
     329        13212 :             ELSE IF (my_neighbor_image) THEN
     330              :                ! Keep the neighbor-list image while selecting one atom-pair representative.
     331           72 :                ra(:) = particle_set(iatom)%r(:)
     332           72 :                rb(:) = ra(:) + rab(:)
     333           18 :                rac(:) = pbc(rc, ra, cell)
     334           72 :                rbc(:) = rac(:) + rab(:)
     335           72 :                dab = NORM2(rab)
     336              :             ELSE
     337              :                ! Fold atomic positions around the reference point.
     338       105552 :                ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
     339       105552 :                rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
     340        52776 :                rab(:) = ra(:) - rb(:)
     341        52776 :                rac(:) = ra(:) - rc(:)
     342        52776 :                rbc(:) = rb(:) - rc(:)
     343        52776 :                dab = NORM2(rab)
     344              :             END IF
     345              : 
     346        70760 :             DO iset = 1, nseta
     347              : 
     348        21054 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
     349        21054 :                sgfa = first_sgfa(1, iset)
     350              : 
     351        73709 :                DO jset = 1, nsetb
     352              : 
     353        39043 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     354              : 
     355        38161 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
     356        38161 :                   sgfb = first_sgfb(1, jset)
     357              : 
     358              :                   ! Calculate the primitive integrals
     359              :                   CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
     360              :                               rpgfa(:, iset), la_min(iset), &
     361              :                               lb_max(jset), npgfb(jset), zetb(:, jset), &
     362        38161 :                               rpgfb(:, jset), nmoments, rac, rbc, mab)
     363              : 
     364              :                   ! Contraction step
     365       179128 :                   DO i = component_start, nm
     366              : 
     367              :                      CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     368              :                                 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
     369              :                                 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     370       119913 :                                 0.0_dp, work(1, 1), SIZE(work, 1))
     371              : 
     372       158956 :                      IF (iatom <= jatom) THEN
     373              : 
     374              :                         CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     375              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     376              :                                    work(1, 1), SIZE(work, 1), &
     377              :                                    1.0_dp, mint(i)%block(sgfa, sgfb), &
     378        78094 :                                    SIZE(mint(i)%block, 1))
     379              : 
     380              :                      ELSE
     381              : 
     382              :                         CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
     383              :                                    1.0_dp, work(1, 1), SIZE(work, 1), &
     384              :                                    sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     385              :                                    1.0_dp, mint(i)%block(sgfb, sgfa), &
     386        41819 :                                    SIZE(mint(i)%block, 1))
     387              : 
     388              :                      END IF
     389              : 
     390              :                   END DO
     391              : 
     392              :                END DO
     393              :             END DO
     394              :          END ASSOCIATE
     395              : 
     396              :       END DO
     397         3428 :       CALL neighbor_list_iterator_release(nl_iterator)
     398              : 
     399              :       ! Release work storage
     400         3428 :       DEALLOCATE (mab, basis_set_list)
     401         3428 :       DEALLOCATE (work)
     402        13882 :       DO i = component_start, nm
     403        13882 :          NULLIFY (mint(i)%block)
     404              :       END DO
     405         3428 :       DEALLOCATE (mint)
     406              : 
     407         3428 :       CALL timestop(handle)
     408              : 
     409         6856 :    END SUBROUTINE build_local_moment_matrix
     410              : 
     411              : ! **************************************************************************************************
     412              : !> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
     413              : !>        Optionally stores the multipole moments themselves for free.
     414              : !>        Note that the multipole moments are symmetric while their derivatives are anti-symmetric
     415              : !>        Only first derivatives are performed, e. g. x d/dy
     416              : !> \param qs_env ...
     417              : !> \param moments_der will contain the derivatives of the multipole moments
     418              : !> \param nmoments_der order of the moments with derivatives
     419              : !> \param nmoments order of the multipole moments (no derivatives, same output as
     420              : !>        build_local_moment_matrix, needs moments as arguments to store results)
     421              : !> \param ref_point ...
     422              : !> \param moments contains the multipole moments, optionally for free, up to order nmoments
     423              : !> \param basis_type basis set used for the integral evaluation
     424              : !> \param minimum_image apply the half-cell image filter
     425              : !> \param ordered use ordered neighbor blocks and, for ordinary derivatives, the full neighbor list
     426              : !> \param lambda optional atom selector for the right-hand derivative
     427              : !> \param deltaR optional weights for left and right nuclear-coordinate derivatives
     428              : !> \param neighbor_image use the neighbor-list image for coordinate construction
     429              : ! **************************************************************************************************
     430         4302 :    SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
     431         4302 :                                              ref_point, moments, basis_type, minimum_image, ordered, lambda, &
     432              :                                              deltaR, neighbor_image)
     433              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     434              :       TYPE(dbcsr_p_type), DIMENSION(:, :), &
     435              :          INTENT(INOUT), POINTER                          :: moments_der
     436              :       INTEGER, INTENT(IN)                                :: nmoments_der, nmoments
     437              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: ref_point
     438              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
     439              :          OPTIONAL, POINTER                               :: moments
     440              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
     441              :       LOGICAL, OPTIONAL                                   :: minimum_image, ordered
     442              :       INTEGER, OPTIONAL                                   :: lambda
     443              :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: deltaR
     444              :       LOGICAL, OPTIONAL                                   :: neighbor_image
     445              : 
     446              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moments_der_matrix'
     447              : 
     448              :       INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
     449              :                  jatom, jkind, jpgf, jset, last_jatom, lda_min, M_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
     450              :                  nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
     451              :       LOGICAL                                            :: accumulate_images, found, my_minimum_image, my_neighbor_image, &
     452              :                                                             my_ordered, new_atom_pair
     453              :       REAL(KIND=dp)                                      :: dab, derivative_sign
     454              :       REAL(KIND=dp), DIMENSION(3)                        :: half_cell
     455         4302 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     456         4302 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: mab
     457         4302 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: difmab
     458              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, rb, rbc, rc
     459         4302 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: mab_tmp
     460         4302 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: mom_block
     461         4302 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :)   :: mom_block_der
     462              :       TYPE(cell_type), POINTER                           :: cell
     463         4302 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     464              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     465              :       TYPE(neighbor_list_iterator_p_type), &
     466         4302 :          DIMENSION(:), POINTER                           :: nl_iterator
     467              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     468         4302 :          POINTER                                         :: sab_all, sab_orb
     469         4302 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     470         4302 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     471              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     472              : 
     473         4302 :       nmom_build = MAX(nmoments, nmoments_der)      ! build moments up to order nmom_buiod
     474         4302 :       IF (nmom_build < 1) RETURN
     475              : 
     476         4302 :       CALL timeset(routineN, handle)
     477              : 
     478         4302 :       nders = 1                                    ! only first order derivatives
     479         4302 :       dimders = ncoset(nders) - 1
     480              : 
     481         4302 :       my_minimum_image = .FALSE.
     482         4302 :       IF (PRESENT(minimum_image)) my_minimum_image = minimum_image
     483         4302 :       my_ordered = .FALSE.
     484         4302 :       IF (PRESENT(ordered)) my_ordered = ordered
     485         4302 :       my_neighbor_image = my_ordered
     486         4302 :       IF (PRESENT(neighbor_image)) my_neighbor_image = neighbor_image
     487         4302 :       CPASSERT(.NOT. (PRESENT(deltaR) .AND. PRESENT(lambda)))
     488         4302 :       IF (PRESENT(lambda) .OR. PRESENT(deltaR)) THEN
     489          510 :          CPASSERT(my_ordered)
     490              :       END IF
     491         4302 :       accumulate_images = my_ordered .AND. my_neighbor_image
     492         4302 :       IF (PRESENT(deltaR) .OR. PRESENT(lambda)) accumulate_images = .FALSE.
     493         4302 :       derivative_sign = 1._dp
     494              :       ! diff_momop returns the right nuclear-center derivative. The former
     495              :       ! rRc_xyz_der_ao interface exposed the corresponding electronic derivative
     496              :       ! on the ordinary ordered path.
     497         4302 :       IF (my_ordered .AND. .NOT. PRESENT(deltaR) .AND. .NOT. PRESENT(lambda)) THEN
     498         3750 :          derivative_sign = -1._dp
     499              :       END IF
     500              : 
     501         4302 :       NULLIFY (qs_kind_set, particle_set, sab_all, sab_orb, cell)
     502              :       CALL get_qs_env(qs_env=qs_env, &
     503              :                       qs_kind_set=qs_kind_set, &
     504              :                       particle_set=particle_set, &
     505              :                       cell=cell, &
     506              :                       sab_orb=sab_orb, &
     507         4302 :                       sab_all=sab_all)
     508              : 
     509         4302 :       IF (my_ordered) THEN
     510         4260 :          CPASSERT(nmoments_der > 0)
     511         4260 :          CPASSERT(dbcsr_get_matrix_type(moments_der(1, 1)%matrix) == dbcsr_type_no_symmetry)
     512         4260 :          IF (nmoments > 0) THEN
     513         3750 :             CPASSERT(PRESENT(moments))
     514         3750 :             CPASSERT(dbcsr_get_matrix_type(moments(1)%matrix) == dbcsr_type_no_symmetry)
     515              :          END IF
     516              :       END IF
     517              : 
     518         4302 :       IF (my_minimum_image) THEN
     519            0 :          half_cell = 0.5_dp*NORM2(cell%hmat, DIM=1)
     520              :       END IF
     521              : 
     522         4302 :       nkind = SIZE(qs_kind_set)
     523              : 
     524              :       ! Work storage
     525              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     526              :                            maxco=maxco, maxsgf=maxsgf, &
     527         4848 :                            basis_type=basis_type)
     528              : 
     529         4302 :       IF (nmoments > 0) THEN
     530         3790 :          nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
     531         3790 :          CPASSERT(SIZE(moments) == nm)
     532              :          ! storage for integrals
     533        18938 :          ALLOCATE (mab(maxco, maxco, nm))
     534              :          ! blocks
     535         3790 :          mab(:, :, :) = 0.0_dp
     536        45480 :          ALLOCATE (mom_block(nm))
     537        37900 :          DO i = 1, nm
     538        37900 :             NULLIFY (mom_block(i)%block)
     539              :          END DO
     540              :       END IF
     541              : 
     542         4302 :       IF (nmoments_der > 0) THEN
     543         4302 :          M_dim = ncoset(nmoments_der) - 1
     544         4302 :          CPASSERT(SIZE(moments_der, dim=1) == M_dim)
     545         4302 :          CPASSERT(SIZE(moments_der, dim=2) == dimders)
     546              :          ! storage for integrals
     547        25800 :          ALLOCATE (difmab(maxco, maxco, M_dim, dimders))
     548         4302 :          difmab(:, :, :, :) = 0.0_dp
     549              :          ! blocks
     550       136368 :          ALLOCATE (mom_block_der(M_dim, dimders))
     551        39720 :          DO i = 1, M_dim
     552       145974 :             DO ider = 1, dimders
     553       141672 :                NULLIFY (mom_block_der(i, ider)%block)
     554              :             END DO
     555              :          END DO
     556              :       END IF
     557              : 
     558        17202 :       ALLOCATE (work(maxco, maxsgf))
     559         4302 :       work(:, :) = 0.0_dp
     560              : 
     561         4302 :       NULLIFY (basis_set_a, basis_set_b, basis_set_list)
     562         4302 :       NULLIFY (qs_kind)
     563        20138 :       ALLOCATE (basis_set_list(nkind))
     564        11534 :       DO ikind = 1, nkind
     565         7232 :          qs_kind => qs_kind_set(ikind)
     566         7232 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
     567        11534 :          IF (ASSOCIATED(basis_set_a)) THEN
     568         7232 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     569              :          ELSE
     570            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     571              :          END IF
     572              :       END DO
     573              : 
     574              :       ! Calculate derivatives looping over neighbour list
     575         4302 :       NULLIFY (nl_iterator)
     576         4302 :       IF (my_ordered) THEN
     577         4260 :          CALL neighbor_list_iterator_create(nl_iterator, sab_all)
     578              :       ELSE
     579           42 :          CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     580              :       END IF
     581       308016 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     582              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
     583       303714 :                                 iatom=iatom, jatom=jatom, r=rab)
     584       303714 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     585       303714 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     586       303714 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     587       303714 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     588              :          ASSOCIATE ( &
     589              :             ! basis ikind
     590              :             first_sgfa => basis_set_a%first_sgf, &
     591              :             la_max => basis_set_a%lmax, &
     592              :             la_min => basis_set_a%lmin, &
     593              :             npgfa => basis_set_a%npgf, &
     594              :             nsgfa => basis_set_a%nsgf_set, &
     595              :             rpgfa => basis_set_a%pgf_radius, &
     596              :             set_radius_a => basis_set_a%set_radius, &
     597              :             sphi_a => basis_set_a%sphi, &
     598              :             zeta => basis_set_a%zet, &
     599              :             ! basis jkind, &
     600              :             first_sgfb => basis_set_b%first_sgf, &
     601              :             lb_max => basis_set_b%lmax, &
     602              :             lb_min => basis_set_b%lmin, &
     603              :             npgfb => basis_set_b%npgf, &
     604              :             nsgfb => basis_set_b%nsgf_set, &
     605              :             rpgfb => basis_set_b%pgf_radius, &
     606              :             set_radius_b => basis_set_b%set_radius, &
     607              :             sphi_b => basis_set_b%sphi, &
     608              :             zetb => basis_set_b%zet)
     609              : 
     610       303714 :             IF (SIZE(sphi_a, 1) == 0) CYCLE
     611       303693 :             IF (SIZE(sphi_b, 1) == 0) CYCLE
     612              : 
     613       303693 :             nseta = basis_set_a%nset
     614       303693 :             nsetb = basis_set_b%nset
     615              : 
     616       303693 :             IF (inode == 1) last_jatom = 0
     617              : 
     618       303693 :             IF (my_minimum_image) THEN
     619            0 :                IF (ANY(ABS(rab(:)) > half_cell(:))) CYCLE
     620              :             END IF
     621              : 
     622       303693 :             new_atom_pair = jatom /= last_jatom
     623              :             ! The ordinary ordered neighbor-image path accumulates all images for one atom pair.
     624              :             ! Specialized derivative paths select one representative image, as in the former builders.
     625       303693 :             IF (.NOT. accumulate_images .AND. .NOT. new_atom_pair) CYCLE
     626              : 
     627       284811 :             IF (new_atom_pair) THEN
     628        21597 :                last_jatom = jatom
     629              : 
     630              :                ! get blocks
     631        21597 :                IF (my_ordered) THEN
     632        20253 :                   irow = iatom
     633        20253 :                   icol = jatom
     634         1344 :                ELSE IF (iatom <= jatom) THEN
     635          710 :                   irow = iatom
     636          710 :                   icol = jatom
     637              :                ELSE
     638          634 :                   irow = jatom
     639          634 :                   icol = iatom
     640              :                END IF
     641              : 
     642        21597 :                IF (nmoments > 0) THEN
     643       192960 :                   DO i = 1, nm
     644       173664 :                      NULLIFY (mom_block(i)%block)
     645              :                      ! get block from pre calculated overlap matrix
     646              :                      CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
     647       173664 :                                             row=irow, col=icol, BLOCK=mom_block(i)%block, found=found)
     648       173664 :                      CPASSERT(found .AND. ASSOCIATED(mom_block(i)%block))
     649     13565520 :                      mom_block(i)%block = 0._dp
     650              :                   END DO
     651              :                END IF
     652        21597 :                IF (nmoments_der > 0) THEN
     653       194172 :                   DO i = 1, M_dim
     654       711897 :                      DO ider = 1, dimders
     655       517725 :                         NULLIFY (mom_block_der(i, ider)%block)
     656              :                         CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
     657              :                                                row=irow, col=icol, &
     658              :                                                block=mom_block_der(i, ider)%block, &
     659       517725 :                                                found=found)
     660       517725 :                         CPASSERT(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
     661     40655619 :                         mom_block_der(i, ider)%block = 0._dp
     662              :                      END DO
     663              :                   END DO
     664              :                END IF
     665              :             END IF
     666              : 
     667              :             ! reference point
     668       284811 :             IF (PRESENT(ref_point)) THEN
     669      1139244 :                rc(:) = ref_point(:)
     670              :             ELSE
     671            0 :                rc(:) = 0._dp
     672              :             END IF
     673              :             ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
     674              :             ! by folding around the center, such screwing can be avoided for a proper choice of center.
     675       284811 :             IF (my_neighbor_image) THEN
     676       281199 :                ra(:) = pbc(particle_set(iatom)%r(:), cell)
     677      1124796 :                rb(:) = ra(:) + rab(:)
     678       281199 :                rac(:) = pbc(rc, ra, cell)
     679      1124796 :                rbc(:) = rac(:) + rab(:)
     680         3612 :             ELSE IF (my_ordered) THEN
     681         9072 :                ra(:) = particle_set(iatom)%r(:)
     682         9072 :                rb(:) = particle_set(jatom)%r(:)
     683         2268 :                rab(:) = pbc(rb, ra, cell)
     684         9072 :                rac(:) = pbc(ra - rc, cell)
     685         9072 :                rbc(:) = pbc(rb - rc, cell)
     686              :             ELSE
     687        10752 :                ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
     688        10752 :                rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
     689              :                ! we dont use PBC at this point
     690         5376 :                rab(:) = ra(:) - rb(:)
     691         5376 :                rac(:) = ra(:) - rc(:)
     692         5376 :                rbc(:) = rb(:) - rc(:)
     693              :             END IF
     694      1139244 :             dab = NORM2(rab)
     695              : 
     696      1243365 :             DO iset = 1, nseta
     697              : 
     698       654840 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
     699       654840 :                sgfa = first_sgfa(1, iset)
     700              : 
     701      2522697 :                DO jset = 1, nsetb
     702              : 
     703      1583046 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     704              : 
     705       600231 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
     706       600231 :                   sgfb = first_sgfb(1, jset)
     707              : 
     708       600231 :                   NULLIFY (mab_tmp)
     709              :                   ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
     710      2974323 :                                     npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
     711              : 
     712              :                   ! Calculate the primitive integrals (need l+1 for derivatives)
     713       600231 :                   lda_min = MAX(0, la_min(iset) - 1)
     714              :                   CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
     715              :                               rpgfa(:, iset), lda_min, &
     716              :                               lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
     717       600231 :                               rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
     718              : 
     719       600231 :                   IF (nmoments_der > 0) THEN
     720              :                      CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
     721              :                                      rpgfa(:, iset), la_min(iset), &
     722              :                                      lb_max(jset), npgfb(jset), zetb(:, jset), &
     723              :                                      rpgfb(:, jset), lb_min(jset), &
     724              :                                      nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp, &
     725       600231 :                                      lambda=lambda, deltaR=deltaR, iatom=iatom, jatom=jatom)
     726              :                   END IF
     727              : 
     728       600231 :                   IF (nmoments > 0) THEN
     729              :                      ! copy subset of mab_tmp (l+1) to mab (l)
     730       597930 :                      mab = 0.0_dp
     731      5979300 :                      DO ii = 1, nm
     732      5381370 :                         na = 0
     733      5381370 :                         nda = 0
     734     17745306 :                         DO ipgf = 1, npgfa(iset)
     735     11766006 :                            nb = 0
     736     11766006 :                            ndb = 0
     737     43098966 :                            DO jpgf = 1, npgfb(jset)
     738    128571651 :                               DO j = 1, ncoset(lb_max(jset))
     739    488676321 :                                  DO i = 1, ncoset(la_max(iset))
     740    457343361 :                                     mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
     741              :                                  END DO ! i
     742              :                               END DO ! j
     743     31332960 :                               nb = nb + ncoset(lb_max(jset))
     744     43098966 :                               ndb = ndb + ncoset(lb_max(jset) + 1)
     745              :                            END DO ! jpgf
     746     11766006 :                            na = na + ncoset(la_max(iset))
     747     17147376 :                            nda = nda + ncoset(la_max(iset) + 1)
     748              :                         END DO ! ipgf
     749              :                      END DO
     750              :                      ! Contraction step
     751      5979300 :                      DO i = 1, nm
     752              : 
     753              :                         CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     754              :                                    1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
     755              :                                    sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     756      5381370 :                                    0.0_dp, work(1, 1), SIZE(work, 1))
     757              : 
     758      5979300 :                         IF (my_ordered .OR. iatom <= jatom) THEN
     759              :                            CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     760              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     761              :                                       work(1, 1), SIZE(work, 1), &
     762              :                                       1.0_dp, mom_block(i)%block(sgfa, sgfb), &
     763      5377653 :                                       SIZE(mom_block(i)%block, 1))
     764              :                         ELSE
     765              :                            CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
     766              :                                       1.0_dp, work(1, 1), SIZE(work, 1), &
     767              :                                       sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     768              :                                       1.0_dp, mom_block(i)%block(sgfb, sgfa), &
     769         3717 :                                       SIZE(mom_block(i)%block, 1))
     770              :                         END IF
     771              :                      END DO
     772              :                   END IF
     773              : 
     774       600231 :                   IF (nmoments_der > 0) THEN
     775      5982852 :                      DO i = 1, M_dim
     776     22130715 :                         DO ider = 1, dimders
     777              :                            CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     778              :                                       derivative_sign, difmab(1, 1, i, ider), SIZE(difmab, 1), &
     779              :                                       sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     780     16147863 :                                       0._dp, work(1, 1), SIZE(work, 1))
     781              : 
     782     21530484 :                            IF (my_ordered .OR. iatom <= jatom) THEN
     783              :                               CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     784              :                                          1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     785              :                                          work(1, 1), SIZE(work, 1), &
     786              :                                          1.0_dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
     787     16144092 :                                          SIZE(mom_block_der(i, ider)%block, 1))
     788              :                            ELSE
     789              :                               CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
     790              :                                          -1.0_dp, work(1, 1), SIZE(work, 1), &
     791              :                                          sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     792              :                                          1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
     793         3771 :                                          SIZE(mom_block_der(i, ider)%block, 1))
     794              :                            END IF
     795              :                         END DO
     796              :                      END DO
     797              :                   END IF
     798      2237886 :                   DEALLOCATE (mab_tmp)
     799              :                END DO
     800              :             END DO
     801              :          END ASSOCIATE
     802              :       END DO
     803         4302 :       CALL neighbor_list_iterator_release(nl_iterator)
     804              : 
     805              :       ! deallocations
     806         4302 :       DEALLOCATE (basis_set_list)
     807         4302 :       DEALLOCATE (work)
     808         4302 :       IF (nmoments > 0) THEN
     809         3790 :          DEALLOCATE (mab)
     810        37900 :          DO i = 1, nm
     811        37900 :             NULLIFY (mom_block(i)%block)
     812              :          END DO
     813         3790 :          DEALLOCATE (mom_block)
     814              :       END IF
     815         4302 :       IF (nmoments_der > 0) THEN
     816         4302 :          DEALLOCATE (difmab)
     817        39720 :          DO i = 1, M_dim
     818       145974 :             DO ider = 1, dimders
     819       141672 :                NULLIFY (mom_block_der(i, ider)%block)
     820              :             END DO
     821              :          END DO
     822         4302 :          DEALLOCATE (mom_block_der)
     823              :       END IF
     824              : 
     825         4302 :       CALL timestop(handle)
     826              : 
     827         8604 :    END SUBROUTINE build_local_moments_der_matrix
     828              : 
     829              : ! **************************************************************************************************
     830              : !> \brief ...
     831              : !> \param qs_env ...
     832              : !> \param magmom ...
     833              : !> \param nmoments ...
     834              : !> \param ref_point ...
     835              : !> \param ref_points ...
     836              : !> \param basis_type ...
     837              : ! **************************************************************************************************
     838           64 :    SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
     839              : 
     840              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     841              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: magmom
     842              :       INTEGER, INTENT(IN)                                :: nmoments
     843              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: ref_point
     844              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     845              :          OPTIONAL                                        :: ref_points
     846              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
     847              : 
     848              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_magmom_matrix'
     849              : 
     850              :       INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
     851              :                  maxsgf, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
     852              :       LOGICAL                                            :: found
     853              :       REAL(KIND=dp)                                      :: dab
     854           64 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     855           64 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: mab
     856              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, rb, rbc, rc
     857           64 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: mint
     858              :       TYPE(cell_type), POINTER                           :: cell
     859           64 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     860              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     861              :       TYPE(neighbor_list_iterator_p_type), &
     862           64 :          DIMENSION(:), POINTER                           :: nl_iterator
     863              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     864           64 :          POINTER                                         :: sab_orb
     865           64 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     866           64 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     867              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     868              : 
     869           64 :       IF (nmoments < 1) RETURN
     870              : 
     871           64 :       CALL timeset(routineN, handle)
     872              : 
     873              :       ! magnetic dipoles/angular moments only
     874           64 :       nm = 3
     875              : 
     876           64 :       NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
     877              :       CALL get_qs_env(qs_env=qs_env, &
     878              :                       qs_kind_set=qs_kind_set, &
     879              :                       particle_set=particle_set, cell=cell, &
     880           64 :                       sab_orb=sab_orb)
     881              : 
     882           64 :       nkind = SIZE(qs_kind_set)
     883              : 
     884              :       ! Allocate work storage
     885              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     886           64 :                            maxco=maxco, maxsgf=maxsgf)
     887              : 
     888          320 :       ALLOCATE (mab(maxco, maxco, nm))
     889           64 :       mab(:, :, :) = 0.0_dp
     890              : 
     891          256 :       ALLOCATE (work(maxco, maxsgf))
     892           64 :       work(:, :) = 0.0_dp
     893              : 
     894          256 :       ALLOCATE (mint(nm))
     895          256 :       DO i = 1, nm
     896          256 :          NULLIFY (mint(i)%block)
     897              :       END DO
     898              : 
     899          356 :       ALLOCATE (basis_set_list(nkind))
     900          228 :       DO ikind = 1, nkind
     901          164 :          qs_kind => qs_kind_set(ikind)
     902          328 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
     903          228 :          IF (ASSOCIATED(basis_set_a)) THEN
     904          164 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     905              :          ELSE
     906            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     907              :          END IF
     908              :       END DO
     909           64 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     910         7696 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     911              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
     912         7632 :                                 iatom=iatom, jatom=jatom, r=rab)
     913         7632 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     914         7632 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     915         7632 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     916         7632 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     917              :          ASSOCIATE ( &
     918              :             ! basis ikind
     919              :             first_sgfa => basis_set_a%first_sgf, &
     920              :             la_max => basis_set_a%lmax, &
     921              :             la_min => basis_set_a%lmin, &
     922              :             npgfa => basis_set_a%npgf, &
     923              :             nsgfa => basis_set_a%nsgf_set, &
     924              :             rpgfa => basis_set_a%pgf_radius, &
     925              :             set_radius_a => basis_set_a%set_radius, &
     926              :             sphi_a => basis_set_a%sphi, &
     927              :             zeta => basis_set_a%zet, &
     928              :             ! basis jkind, &
     929              :             first_sgfb => basis_set_b%first_sgf, &
     930              :             lb_max => basis_set_b%lmax, &
     931              :             lb_min => basis_set_b%lmin, &
     932              :             npgfb => basis_set_b%npgf, &
     933              :             nsgfb => basis_set_b%nsgf_set, &
     934              :             rpgfb => basis_set_b%pgf_radius, &
     935              :             set_radius_b => basis_set_b%set_radius, &
     936              :             sphi_b => basis_set_b%sphi, &
     937              :             zetb => basis_set_b%zet)
     938              : 
     939         7632 :             nseta = basis_set_a%nset
     940         7632 :             nsetb = basis_set_b%nset
     941              : 
     942         7632 :             IF (iatom <= jatom) THEN
     943         4216 :                irow = iatom
     944         4216 :                icol = jatom
     945              :             ELSE
     946         3416 :                irow = jatom
     947         3416 :                icol = iatom
     948              :             END IF
     949              : 
     950        30528 :             DO i = 1, nm
     951        22896 :                NULLIFY (mint(i)%block)
     952              :                CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
     953        22896 :                                       row=irow, col=icol, BLOCK=mint(i)%block, found=found)
     954        22896 :                CPASSERT(found .AND. ASSOCIATED(mint(i)%block))
     955       384762 :                mint(i)%block = 0._dp
     956              :             END DO
     957              : 
     958              :             ! fold atomic position back into unit cell
     959         7632 :             IF (PRESENT(ref_points)) THEN
     960            0 :                rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
     961         7632 :             ELSE IF (PRESENT(ref_point)) THEN
     962        30528 :                rc(:) = ref_point(:)
     963              :             ELSE
     964            0 :                rc(:) = 0._dp
     965              :             END IF
     966              :             ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
     967              :             ! by folding around the center, such screwing can be avoided for a proper choice of center.
     968        61056 :             ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
     969        61056 :             rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
     970              :             ! we dont use PBC at this point
     971        30528 :             rab(:) = ra(:) - rb(:)
     972        30528 :             rac(:) = ra(:) - rc(:)
     973        30528 :             rbc(:) = rb(:) - rc(:)
     974        30528 :             dab = NORM2(rab)
     975              : 
     976        23013 :             DO iset = 1, nseta
     977              : 
     978         7749 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
     979         7749 :                sgfa = first_sgfa(1, iset)
     980              : 
     981        23364 :                DO jset = 1, nsetb
     982              : 
     983         7983 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     984              : 
     985         6990 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
     986         6990 :                   sgfb = first_sgfb(1, jset)
     987              : 
     988              :                   ! Calculate the primitive integrals
     989              :                   CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
     990              :                               rpgfa(:, iset), la_min(iset), &
     991              :                               lb_max(jset), npgfb(jset), zetb(:, jset), &
     992         6990 :                               rpgfb(:, jset), rac, rbc, mab)
     993              : 
     994              :                   ! Contraction step
     995        35709 :                   DO i = 1, nm
     996              :                      CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     997              :                                 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
     998              :                                 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     999        20970 :                                 0.0_dp, work(1, 1), SIZE(work, 1))
    1000              : 
    1001        28953 :                      IF (iatom <= jatom) THEN
    1002              :                         CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
    1003              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1004              :                                    work(1, 1), SIZE(work, 1), &
    1005              :                                    1.0_dp, mint(i)%block(sgfa, sgfb), &
    1006        11901 :                                    SIZE(mint(i)%block, 1))
    1007              :                      ELSE
    1008              :                         CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
    1009              :                                    -1.0_dp, work(1, 1), SIZE(work, 1), &
    1010              :                                    sphi_a(1, sgfa), SIZE(sphi_a, 1), &
    1011              :                                    1.0_dp, mint(i)%block(sgfb, sgfa), &
    1012         9069 :                                    SIZE(mint(i)%block, 1))
    1013              :                      END IF
    1014              : 
    1015              :                   END DO
    1016              : 
    1017              :                END DO
    1018              :             END DO
    1019              :          END ASSOCIATE
    1020              :       END DO
    1021           64 :       CALL neighbor_list_iterator_release(nl_iterator)
    1022              : 
    1023              :       ! Release work storage
    1024           64 :       DEALLOCATE (mab, basis_set_list)
    1025           64 :       DEALLOCATE (work)
    1026          256 :       DO i = 1, nm
    1027          256 :          NULLIFY (mint(i)%block)
    1028              :       END DO
    1029           64 :       DEALLOCATE (mint)
    1030              : 
    1031           64 :       CALL timestop(handle)
    1032              : 
    1033          192 :    END SUBROUTINE build_local_magmom_matrix
    1034              : 
    1035              : ! **************************************************************************************************
    1036              : !> \brief ...
    1037              : !> \param qs_env ...
    1038              : !> \param cosmat ...
    1039              : !> \param sinmat ...
    1040              : !> \param kvec ...
    1041              : !> \param sab_orb_external ...
    1042              : !> \param basis_type ...
    1043              : ! **************************************************************************************************
    1044        13456 :    SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
    1045              : 
    1046              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1047              :       TYPE(dbcsr_type), POINTER                          :: cosmat, sinmat
    1048              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: kvec
    1049              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1050              :          OPTIONAL, POINTER                               :: sab_orb_external
    1051              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
    1052              : 
    1053              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_moment_matrix'
    1054              : 
    1055              :       INTEGER                                            :: handle
    1056        53824 :       TYPE(dbcsr_p_type), DIMENSION(2, 1)                :: op_sm_set
    1057              :       REAL(KIND=dp), DIMENSION(3, 1)                     :: kvec_batch
    1058              : 
    1059        13456 :       CALL timeset(routineN, handle)
    1060              : 
    1061        13456 :       op_sm_set(1, 1)%matrix => cosmat
    1062        13456 :       op_sm_set(2, 1)%matrix => sinmat
    1063        13456 :       kvec_batch(:, 1) = kvec(:)
    1064              : 
    1065              :       CALL build_exp_ikr_matrix(qs_env, op_sm_set, kvec_batch, &
    1066        15752 :                                 sab_orb_external=sab_orb_external, basis_type=basis_type)
    1067              : 
    1068        13456 :       CALL timestop(handle)
    1069              : 
    1070        13456 :    END SUBROUTINE build_berry_moment_matrix
    1071              : 
    1072              : ! **************************************************************************************************
    1073              : !> \brief ...
    1074              : !> \param qs_env ...
    1075              : !> \param cosmat ...
    1076              : !> \param sinmat ...
    1077              : !> \param kvec ...
    1078              : !> \param basis_type ...
    1079              : !> \param ordered use all ordered atom pairs for a non-Hermitian cross-k overlap
    1080              : !> \param ordered_kpoints image-cell mapping for the ordered neighbour list
    1081              : ! **************************************************************************************************
    1082          116 :    SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type, ordered, ordered_kpoints)
    1083              : 
    1084              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1085              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: cosmat, sinmat
    1086              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: kvec
    1087              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
    1088              :       LOGICAL, INTENT(IN), OPTIONAL                      :: ordered
    1089              :       TYPE(kpoint_type), POINTER, OPTIONAL               :: ordered_kpoints
    1090              : 
    1091              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_kpoint_matrix'
    1092              : 
    1093              :       CHARACTER                                          :: matrix_symmetry
    1094              :       INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
    1095              :                  ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
    1096              :       INTEGER, DIMENSION(3)                              :: icell
    1097          116 :       INTEGER, DIMENSION(:), POINTER                     :: row_blk_sizes
    1098          116 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    1099              :       LOGICAL                                            :: found, full_pairs, use_cell_mapping
    1100          116 :       REAL(dp), DIMENSION(:, :), POINTER                 :: cblock, cosab, sblock, sinab, work
    1101              :       REAL(KIND=dp)                                      :: dab
    1102              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rb
    1103              :       TYPE(cell_type), POINTER                           :: cell
    1104              :       TYPE(dbcsr_distribution_type), POINTER             :: dbcsr_dist
    1105              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1106          116 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
    1107              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set, basis_set_a, basis_set_b
    1108              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1109              :       TYPE(neighbor_list_iterator_p_type), &
    1110          116 :          DIMENSION(:), POINTER                           :: nl_iterator
    1111              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1112          116 :          POINTER                                         :: sab_orb
    1113          116 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1114          116 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1115              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    1116              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1117              : 
    1118          116 :       CALL timeset(routineN, handle)
    1119              : 
    1120          116 :       full_pairs = .FALSE.
    1121          116 :       IF (PRESENT(ordered)) full_pairs = ordered
    1122          116 :       matrix_symmetry = dbcsr_type_symmetric
    1123          116 :       IF (full_pairs) matrix_symmetry = dbcsr_type_no_symmetry
    1124              : 
    1125              :       CALL get_qs_env(qs_env, &
    1126              :                       ks_env=ks_env, &
    1127          116 :                       dft_control=dft_control)
    1128          116 :       nimg = dft_control%nimages
    1129          116 :       IF (nimg > 1) THEN
    1130          116 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
    1131          116 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
    1132          116 :          use_cell_mapping = .TRUE.
    1133              :       ELSE
    1134              :          use_cell_mapping = .FALSE.
    1135              :       END IF
    1136              : 
    1137              :       CALL get_qs_env(qs_env=qs_env, &
    1138              :                       qs_kind_set=qs_kind_set, &
    1139              :                       particle_set=particle_set, cell=cell, &
    1140          116 :                       sab_orb=sab_orb)
    1141              : 
    1142          116 :       IF (full_pairs) THEN
    1143           20 :          CALL get_qs_env(qs_env, sab_all=sab_orb)
    1144           20 :          IF (.NOT. ASSOCIATED(sab_orb)) CPABORT("Ordered Berry integrals require k-point neighbour lists.")
    1145           20 :          CPASSERT(PRESENT(ordered_kpoints))
    1146           20 :          CALL get_kpoint_info(ordered_kpoints, cell_to_index=cell_to_index)
    1147           20 :          nimg = SIZE(ordered_kpoints%index_to_cell, 2)
    1148           20 :          use_cell_mapping = .TRUE.
    1149              :       END IF
    1150              : 
    1151          116 :       nkind = SIZE(qs_kind_set)
    1152          116 :       natom = SIZE(particle_set)
    1153          464 :       ALLOCATE (basis_set_list(nkind))
    1154          232 :       DO ikind = 1, nkind
    1155          116 :          qs_kind => qs_kind_set(ikind)
    1156          232 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
    1157          232 :          IF (ASSOCIATED(basis_set)) THEN
    1158          116 :             basis_set_list(ikind)%gto_basis_set => basis_set
    1159              :          ELSE
    1160            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
    1161              :          END IF
    1162              :       END DO
    1163              : 
    1164          348 :       ALLOCATE (row_blk_sizes(natom))
    1165              :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
    1166          116 :                             basis=basis_set_list)
    1167          116 :       CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
    1168              :       ! (re)allocate matrix sets
    1169          116 :       CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
    1170          116 :       CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
    1171        17582 :       DO i = 1, nimg
    1172              :          ! sin
    1173        17466 :          ALLOCATE (sinmat(1, i)%matrix)
    1174              :          CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
    1175              :                            name="SINMAT", &
    1176              :                            dist=dbcsr_dist, matrix_type=matrix_symmetry, &
    1177        17466 :                            row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
    1178        17466 :          CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
    1179        17466 :          CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
    1180              :          ! cos
    1181        17466 :          ALLOCATE (cosmat(1, i)%matrix)
    1182              :          CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
    1183              :                            name="COSMAT", &
    1184              :                            dist=dbcsr_dist, matrix_type=matrix_symmetry, &
    1185        17466 :                            row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
    1186        17466 :          CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
    1187        17582 :          CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
    1188              :       END DO
    1189              : 
    1190          116 :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
    1191          116 :       ldab = ldwork
    1192          464 :       ALLOCATE (cosab(ldab, ldab))
    1193          348 :       ALLOCATE (sinab(ldab, ldab))
    1194          348 :       ALLOCATE (work(ldwork, ldwork))
    1195              : 
    1196          116 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
    1197        98908 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1198              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
    1199        98792 :                                 iatom=iatom, jatom=jatom, r=rab, cell=icell)
    1200        98792 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
    1201        98792 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    1202        98792 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
    1203        98792 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    1204              :          ASSOCIATE ( &
    1205              :             ! basis ikind
    1206              :             first_sgfa => basis_set_a%first_sgf, &
    1207              :             la_max => basis_set_a%lmax, &
    1208              :             la_min => basis_set_a%lmin, &
    1209              :             npgfa => basis_set_a%npgf, &
    1210              :             nsgfa => basis_set_a%nsgf_set, &
    1211              :             rpgfa => basis_set_a%pgf_radius, &
    1212              :             set_radius_a => basis_set_a%set_radius, &
    1213              :             sphi_a => basis_set_a%sphi, &
    1214              :             zeta => basis_set_a%zet, &
    1215              :             ! basis jkind, &
    1216              :             first_sgfb => basis_set_b%first_sgf, &
    1217              :             lb_max => basis_set_b%lmax, &
    1218              :             lb_min => basis_set_b%lmin, &
    1219              :             npgfb => basis_set_b%npgf, &
    1220              :             nsgfb => basis_set_b%nsgf_set, &
    1221              :             rpgfb => basis_set_b%pgf_radius, &
    1222              :             set_radius_b => basis_set_b%set_radius, &
    1223              :             sphi_b => basis_set_b%sphi, &
    1224       197584 :             zetb => basis_set_b%zet)
    1225              : 
    1226        98792 :             nseta = basis_set_a%nset
    1227        98792 :             nsetb = basis_set_b%nset
    1228              : 
    1229        98792 :             ldsa = SIZE(sphi_a, 1)
    1230        98792 :             ldsb = SIZE(sphi_b, 1)
    1231              : 
    1232        42450 :             IF (iatom <= jatom .OR. full_pairs) THEN
    1233        56552 :                irow = iatom
    1234        56552 :                icol = jatom
    1235              :             ELSE
    1236        42240 :                irow = jatom
    1237        42240 :                icol = iatom
    1238              :             END IF
    1239              : 
    1240        98792 :             IF (use_cell_mapping) THEN
    1241        98792 :                ic = cell_to_index(icell(1), icell(2), icell(3))
    1242        98792 :                CPASSERT(ic > 0)
    1243              :             ELSE
    1244              :                ic = 1
    1245              :             END IF
    1246              : 
    1247        98792 :             NULLIFY (sblock)
    1248              :             CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
    1249        98792 :                                    row=irow, col=icol, BLOCK=sblock, found=found)
    1250        98792 :             CPASSERT(found)
    1251        98792 :             NULLIFY (cblock)
    1252              :             CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
    1253        98792 :                                    row=irow, col=icol, BLOCK=cblock, found=found)
    1254        98792 :             CPASSERT(found)
    1255              : 
    1256        98792 :             ra(:) = pbc(particle_set(iatom)%r(:), cell)
    1257       395168 :             rb(:) = ra + rab
    1258        98792 :             dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
    1259              : 
    1260       300038 :             DO iset = 1, nseta
    1261              : 
    1262       102454 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
    1263       102454 :                sgfa = first_sgfa(1, iset)
    1264              : 
    1265       311024 :                DO jset = 1, nsetb
    1266              : 
    1267       109778 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
    1268              : 
    1269       106986 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
    1270       106986 :                   sgfb = first_sgfb(1, jset)
    1271              : 
    1272              :                   ! Calculate the primitive integrals
    1273              :                   CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
    1274              :                               lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
    1275       106986 :                               ra, rb, kvec, cosab, sinab)
    1276              :                   CALL contract_cossin(cblock, sblock, &
    1277              :                                        iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
    1278              :                                        jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
    1279       212232 :                                        cosab, sinab, ldab, work, ldwork, ordered=full_pairs)
    1280              : 
    1281              :                END DO
    1282              :             END DO
    1283              :          END ASSOCIATE
    1284              :       END DO
    1285          116 :       CALL neighbor_list_iterator_release(nl_iterator)
    1286              : 
    1287          116 :       DEALLOCATE (cosab)
    1288          116 :       DEALLOCATE (sinab)
    1289          116 :       DEALLOCATE (work)
    1290          116 :       DEALLOCATE (basis_set_list)
    1291          116 :       DEALLOCATE (row_blk_sizes)
    1292              : 
    1293          116 :       CALL timestop(handle)
    1294              : 
    1295          232 :    END SUBROUTINE build_berry_kpoint_matrix
    1296              : 
    1297              : ! **************************************************************************************************
    1298              : !> \brief ...
    1299              : !> \param qs_env ...
    1300              : !> \param magnetic ...
    1301              : !> \param nmoments ...
    1302              : !> \param reference ...
    1303              : !> \param ref_point ...
    1304              : !> \param unit_number ...
    1305              : ! **************************************************************************************************
    1306          478 :    SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
    1307              : 
    1308              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1309              :       LOGICAL, INTENT(IN)                                :: magnetic
    1310              :       INTEGER, INTENT(IN)                                :: nmoments, reference
    1311              :       REAL(dp), DIMENSION(:), POINTER                    :: ref_point
    1312              :       INTEGER, INTENT(IN)                                :: unit_number
    1313              : 
    1314              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_moment_berry_phase'
    1315              : 
    1316          478 :       CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:)        :: rlab
    1317              :       CHARACTER(LEN=default_string_length)               :: description
    1318              :       COMPLEX(dp)                                        :: xphase(3), zdet, zdeta, zi(3), &
    1319              :                                                             zij(3, 3), zijk(3, 3, 3), &
    1320              :                                                             zijkl(3, 3, 3, 3), zphase(3), zz
    1321              :       INTEGER                                            :: handle, i, ia, idim, ikind, ispin, ix, &
    1322              :                                                             iy, iz, j, k, l, nao, nm, nmo, nmom, &
    1323              :                                                             nmotot, tmp_dim
    1324              :       LOGICAL                                            :: floating, ghost, uniform
    1325              :       REAL(dp)                                           :: charge, ci(3), cij(3, 3), dd, occ, trace
    1326          478 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: mmom
    1327          478 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: rmom
    1328              :       REAL(dp), DIMENSION(3)                             :: kvec, qq, rcc, ria
    1329              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1330              :       TYPE(cell_type), POINTER                           :: cell
    1331          478 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: eigrmat
    1332              :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
    1333          478 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: opvec
    1334          478 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: op_fm_set
    1335          478 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
    1336              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1337              :       TYPE(cp_result_type), POINTER                      :: results
    1338          478 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao
    1339              :       TYPE(dbcsr_type), POINTER                          :: cosmat, sinmat
    1340              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1341              :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1342          478 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1343              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1344          478 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1345          478 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1346              :       TYPE(qs_rho_type), POINTER                         :: rho
    1347              :       TYPE(rt_prop_type), POINTER                        :: rtp
    1348              : 
    1349            0 :       CPASSERT(ASSOCIATED(qs_env))
    1350              : 
    1351          478 :       IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
    1352            0 :          IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
    1353            0 :          RETURN
    1354              :       END IF
    1355              : 
    1356          478 :       CALL timeset(routineN, handle)
    1357              : 
    1358              :       ! restrict maximum moment available
    1359          478 :       nmom = MIN(nmoments, 2)
    1360              : 
    1361          478 :       nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
    1362              :       ! rmom(:,1)=electronic
    1363              :       ! rmom(:,2)=nuclear
    1364              :       ! rmom(:,1)=total
    1365         2390 :       ALLOCATE (rmom(nm + 1, 3))
    1366         1434 :       ALLOCATE (rlab(nm + 1))
    1367          478 :       rmom = 0.0_dp
    1368         2390 :       rlab = ""
    1369          478 :       IF (magnetic) THEN
    1370            0 :          nm = 3
    1371            0 :          ALLOCATE (mmom(nm))
    1372            0 :          mmom = 0._dp
    1373              :       END IF
    1374              : 
    1375          478 :       NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
    1376          478 :                local_particles, matrix_s, mos, rho_ao)
    1377              : 
    1378              :       CALL get_qs_env(qs_env, &
    1379              :                       dft_control=dft_control, &
    1380              :                       rho=rho, &
    1381              :                       cell=cell, &
    1382              :                       results=results, &
    1383              :                       particle_set=particle_set, &
    1384              :                       qs_kind_set=qs_kind_set, &
    1385              :                       para_env=para_env, &
    1386              :                       local_particles=local_particles, &
    1387              :                       matrix_s=matrix_s, &
    1388          478 :                       mos=mos)
    1389              : 
    1390          478 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
    1391              : 
    1392          478 :       NULLIFY (cosmat, sinmat)
    1393          478 :       ALLOCATE (cosmat, sinmat)
    1394          478 :       CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
    1395          478 :       CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
    1396              : 
    1397         2934 :       ALLOCATE (op_fm_set(2, dft_control%nspins))
    1398         1934 :       ALLOCATE (opvec(dft_control%nspins))
    1399         1934 :       ALLOCATE (eigrmat(dft_control%nspins))
    1400          478 :       nmotot = 0
    1401          978 :       DO ispin = 1, dft_control%nspins
    1402          500 :          NULLIFY (tmp_fm_struct, mo_coeff)
    1403          500 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
    1404          500 :          nmotot = nmotot + nmo
    1405          500 :          CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
    1406              :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
    1407          500 :                                   ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
    1408         1500 :          DO i = 1, SIZE(op_fm_set, 1)
    1409         1500 :             CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
    1410              :          END DO
    1411          500 :          CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
    1412         1478 :          CALL cp_fm_struct_release(tmp_fm_struct)
    1413              :       END DO
    1414              : 
    1415              :       ! occupation
    1416          978 :       DO ispin = 1, dft_control%nspins
    1417          500 :          CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
    1418          978 :          IF (.NOT. uniform) THEN
    1419            0 :             CPWARN("Berry phase moments for non uniform MOs' occupation numbers not implemented")
    1420              :          END IF
    1421              :       END DO
    1422              : 
    1423              :       ! reference point
    1424          478 :       CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
    1425         1912 :       rcc = pbc(rcc, cell)
    1426              : 
    1427              :       ! label
    1428         1912 :       DO l = 1, nm
    1429         1434 :          ix = indco(1, l + 1)
    1430         1434 :          iy = indco(2, l + 1)
    1431         1434 :          iz = indco(3, l + 1)
    1432         1912 :          CALL set_label(rlab(l + 1), ix, iy, iz)
    1433              :       END DO
    1434              : 
    1435              :       ! nuclear contribution
    1436         1932 :       DO ia = 1, SIZE(particle_set)
    1437         1454 :          atomic_kind => particle_set(ia)%atomic_kind
    1438         1454 :          CALL get_atomic_kind(atomic_kind, kind_number=ikind)
    1439         1454 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
    1440         1932 :          IF (.NOT. ghost .AND. .NOT. floating) THEN
    1441         1454 :             rmom(1, 2) = rmom(1, 2) - charge
    1442              :          END IF
    1443              :       END DO
    1444         7648 :       ria = twopi*MATMUL(cell%h_inv, rcc)
    1445         1912 :       zphase = CMPLX(COS(ria), SIN(ria), dp)**rmom(1, 2)
    1446              : 
    1447          478 :       zi = 0._dp
    1448          478 :       zij = 0._dp
    1449              :       zijk = 0._dp
    1450              :       zijkl = 0._dp
    1451              : 
    1452          956 :       DO l = 1, nmom
    1453          478 :          SELECT CASE (l)
    1454              :          CASE (1)
    1455              :             ! Dipole
    1456         1912 :             zi(:) = CMPLX(1._dp, 0._dp, dp)
    1457         1932 :             DO ia = 1, SIZE(particle_set)
    1458         1454 :                atomic_kind => particle_set(ia)%atomic_kind
    1459         1454 :                CALL get_atomic_kind(atomic_kind, kind_number=ikind)
    1460         1454 :                CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
    1461         1932 :                IF (.NOT. ghost .AND. .NOT. floating) THEN
    1462         5816 :                   ria = particle_set(ia)%r
    1463         5816 :                   ria = pbc(ria, cell)
    1464         5816 :                   DO i = 1, 3
    1465        17448 :                      kvec(:) = twopi*cell%h_inv(i, :)
    1466        17448 :                      dd = SUM(kvec(:)*ria(:))
    1467         4362 :                      zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
    1468         5816 :                      zi(i) = zi(i)*zdeta
    1469              :                   END DO
    1470              :                END IF
    1471              :             END DO
    1472         1912 :             zi = zi*zphase
    1473         1912 :             ci = AIMAG(LOG(zi))/twopi
    1474         1912 :             qq = AIMAG(LOG(zi))
    1475         7648 :             rmom(2:4, 2) = MATMUL(cell%hmat, ci)
    1476              :          CASE (2)
    1477              :             ! Quadrupole
    1478            0 :             CPABORT("Berry phase moments bigger than 1 not implemented")
    1479            0 :             zij(:, :) = CMPLX(1._dp, 0._dp, dp)
    1480            0 :             DO ia = 1, SIZE(particle_set)
    1481            0 :                atomic_kind => particle_set(ia)%atomic_kind
    1482            0 :                CALL get_atomic_kind(atomic_kind, kind_number=ikind)
    1483            0 :                CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
    1484            0 :                ria = particle_set(ia)%r
    1485            0 :                ria = pbc(ria, cell)
    1486            0 :                DO i = 1, 3
    1487            0 :                   DO j = i, 3
    1488            0 :                      kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
    1489            0 :                      dd = SUM(kvec(:)*ria(:))
    1490            0 :                      zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge
    1491            0 :                      zij(i, j) = zij(i, j)*zdeta
    1492            0 :                      zij(j, i) = zij(i, j)
    1493              :                   END DO
    1494              :                END DO
    1495              :             END DO
    1496            0 :             DO i = 1, 3
    1497            0 :                DO j = 1, 3
    1498            0 :                   zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
    1499            0 :                   zz = zij(i, j)/zi(i)/zi(j)
    1500            0 :                   cij(i, j) = AIMAG(LOG(zz))/twopi
    1501              :                END DO
    1502              :             END DO
    1503            0 :             cij = 0.5_dp*cij/twopi/twopi
    1504            0 :             cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
    1505            0 :             DO k = 4, 9
    1506            0 :                ix = indco(1, k + 1)
    1507            0 :                iy = indco(2, k + 1)
    1508            0 :                iz = indco(3, k + 1)
    1509            0 :                IF (ix == 0) THEN
    1510            0 :                   rmom(k + 1, 2) = cij(iy, iz)
    1511            0 :                ELSE IF (iy == 0) THEN
    1512            0 :                   rmom(k + 1, 2) = cij(ix, iz)
    1513            0 :                ELSE IF (iz == 0) THEN
    1514            0 :                   rmom(k + 1, 2) = cij(ix, iy)
    1515              :                END IF
    1516              :             END DO
    1517              :          CASE (3)
    1518              :             ! Octapole
    1519            0 :             CPABORT("Berry phase moments bigger than 2 not implemented")
    1520              :          CASE (4)
    1521              :             ! Hexadecapole
    1522            0 :             CPABORT("Berry phase moments bigger than 3 not implemented")
    1523              :          CASE DEFAULT
    1524          478 :             CPABORT("Berry phase moments bigger than 4 not implemented")
    1525              :          END SELECT
    1526              :       END DO
    1527              : 
    1528              :       ! electronic contribution
    1529              : 
    1530         7648 :       ria = twopi*REAL(nmotot, dp)*occ*MATMUL(cell%h_inv, rcc)
    1531         1912 :       xphase = CMPLX(COS(ria), SIN(ria), dp)
    1532              : 
    1533              :       ! charge
    1534          478 :       trace = 0.0_dp
    1535          978 :       DO ispin = 1, dft_control%nspins
    1536          500 :          CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
    1537          978 :          rmom(1, 1) = rmom(1, 1) + trace
    1538              :       END DO
    1539              : 
    1540          478 :       zi = 0._dp
    1541          478 :       zij = 0._dp
    1542              :       zijk = 0._dp
    1543              :       zijkl = 0._dp
    1544              : 
    1545          956 :       DO l = 1, nmom
    1546          478 :          SELECT CASE (l)
    1547              :          CASE (1)
    1548              :             ! Dipole
    1549         1912 :             DO i = 1, 3
    1550         5736 :                kvec(:) = twopi*cell%h_inv(i, :)
    1551         1434 :                CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
    1552         1434 :                IF (qs_env%run_rtp) THEN
    1553           48 :                   CALL get_qs_env(qs_env, rtp=rtp)
    1554           48 :                   CALL get_rtp(rtp, mos_new=mos_new)
    1555           48 :                   CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
    1556              :                ELSE
    1557         1386 :                   CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
    1558              :                END IF
    1559         1434 :                zdet = CMPLX(1._dp, 0._dp, dp)
    1560         2934 :                DO ispin = 1, dft_control%nspins
    1561         1500 :                   CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
    1562         8322 :                   DO idim = 1, tmp_dim
    1563              :                      eigrmat(ispin)%local_data(:, idim) = &
    1564              :                         CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
    1565        43968 :                               -op_fm_set(2, ispin)%local_data(:, idim), dp)
    1566              :                   END DO
    1567              :                   ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
    1568         1500 :                   CALL cp_cfm_det(eigrmat(ispin), zdeta)
    1569         1500 :                   zdet = zdet*zdeta
    1570         4434 :                   IF (dft_control%nspins == 1) THEN
    1571         1368 :                      zdet = zdet*zdeta
    1572              :                   END IF
    1573              :                END DO
    1574         1912 :                zi(i) = zdet
    1575              :             END DO
    1576         1912 :             zi = zi*xphase
    1577              :          CASE (2)
    1578              :             ! Quadrupole
    1579            0 :             CPABORT("Berry phase moments bigger than 1 not implemented")
    1580            0 :             DO i = 1, 3
    1581            0 :                DO j = i, 3
    1582            0 :                   kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
    1583            0 :                   CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
    1584            0 :                   IF (qs_env%run_rtp) THEN
    1585            0 :                      CALL get_qs_env(qs_env, rtp=rtp)
    1586            0 :                      CALL get_rtp(rtp, mos_new=mos_new)
    1587            0 :                      CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
    1588              :                   ELSE
    1589            0 :                      CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
    1590              :                   END IF
    1591            0 :                   zdet = CMPLX(1._dp, 0._dp, dp)
    1592            0 :                   DO ispin = 1, dft_control%nspins
    1593            0 :                      CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
    1594            0 :                      DO idim = 1, tmp_dim
    1595              :                         eigrmat(ispin)%local_data(:, idim) = &
    1596              :                            CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
    1597            0 :                                  -op_fm_set(2, ispin)%local_data(:, idim), dp)
    1598              :                      END DO
    1599              :                      ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
    1600            0 :                      CALL cp_cfm_det(eigrmat(ispin), zdeta)
    1601            0 :                      zdet = zdet*zdeta
    1602            0 :                      IF (dft_control%nspins == 1) THEN
    1603            0 :                         zdet = zdet*zdeta
    1604              :                      END IF
    1605              :                   END DO
    1606            0 :                   zij(i, j) = zdet*xphase(i)*xphase(j)
    1607            0 :                   zij(j, i) = zdet*xphase(i)*xphase(j)
    1608              :                END DO
    1609              :             END DO
    1610              :          CASE (3)
    1611              :             ! Octapole
    1612            0 :             CPABORT("Berry phase moments bigger than 2 not implemented")
    1613              :          CASE (4)
    1614              :             ! Hexadecapole
    1615            0 :             CPABORT("Berry phase moments bigger than 3 not implemented")
    1616              :          CASE DEFAULT
    1617          478 :             CPABORT("Berry phase moments bigger than 4 not implemented")
    1618              :          END SELECT
    1619              :       END DO
    1620          956 :       DO l = 1, nmom
    1621          478 :          SELECT CASE (l)
    1622              :          CASE (1)
    1623              :             ! Dipole (apply periodic (2 Pi) boundary conditions)
    1624         1912 :             ci = AIMAG(LOG(zi))
    1625         1912 :             DO i = 1, 3
    1626         1434 :                IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
    1627         1912 :                IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
    1628              :             END DO
    1629         9082 :             rmom(2:4, 1) = MATMUL(cell%hmat, ci)/twopi
    1630              :          CASE (2)
    1631              :             ! Quadrupole
    1632            0 :             CPABORT("Berry phase moments bigger than 1 not implemented")
    1633            0 :             DO i = 1, 3
    1634            0 :                DO j = 1, 3
    1635            0 :                   zz = zij(i, j)/zi(i)/zi(j)
    1636            0 :                   cij(i, j) = AIMAG(LOG(zz))/twopi
    1637              :                END DO
    1638              :             END DO
    1639            0 :             cij = 0.5_dp*cij/twopi/twopi
    1640            0 :             cij = MATMUL(MATMUL(cell%hmat, cij), TRANSPOSE(cell%hmat))
    1641            0 :             DO k = 4, 9
    1642            0 :                ix = indco(1, k + 1)
    1643            0 :                iy = indco(2, k + 1)
    1644            0 :                iz = indco(3, k + 1)
    1645            0 :                IF (ix == 0) THEN
    1646            0 :                   rmom(k + 1, 1) = cij(iy, iz)
    1647            0 :                ELSE IF (iy == 0) THEN
    1648            0 :                   rmom(k + 1, 1) = cij(ix, iz)
    1649            0 :                ELSE IF (iz == 0) THEN
    1650            0 :                   rmom(k + 1, 1) = cij(ix, iy)
    1651              :                END IF
    1652              :             END DO
    1653              :          CASE (3)
    1654              :             ! Octapole
    1655            0 :             CPABORT("Berry phase moments bigger than 2 not implemented")
    1656              :          CASE (4)
    1657              :             ! Hexadecapole
    1658            0 :             CPABORT("Berry phase moments bigger than 3 not implemented")
    1659              :          CASE DEFAULT
    1660          478 :             CPABORT("Berry phase moments bigger than 4 not implemented")
    1661              :          END SELECT
    1662              :       END DO
    1663              : 
    1664         2390 :       rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
    1665          478 :       description = "[DIPOLE]"
    1666          478 :       CALL cp_results_erase(results=results, description=description)
    1667              :       CALL put_results(results=results, description=description, &
    1668          478 :                        values=rmom(2:4, 3))
    1669          478 :       IF (magnetic) THEN
    1670            0 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE., mmom=mmom)
    1671              :       ELSE
    1672          478 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.TRUE.)
    1673              :       END IF
    1674              : 
    1675          478 :       DEALLOCATE (rmom)
    1676          478 :       DEALLOCATE (rlab)
    1677          478 :       IF (magnetic) THEN
    1678            0 :          DEALLOCATE (mmom)
    1679              :       END IF
    1680              : 
    1681          478 :       CALL dbcsr_deallocate_matrix(cosmat)
    1682          478 :       CALL dbcsr_deallocate_matrix(sinmat)
    1683              : 
    1684          478 :       CALL cp_fm_release(opvec)
    1685          478 :       CALL cp_fm_release(op_fm_set)
    1686          978 :       DO ispin = 1, dft_control%nspins
    1687          978 :          CALL cp_cfm_release(eigrmat(ispin))
    1688              :       END DO
    1689          478 :       DEALLOCATE (eigrmat)
    1690              : 
    1691          478 :       CALL timestop(handle)
    1692              : 
    1693         1434 :    END SUBROUTINE qs_moment_berry_phase
    1694              : 
    1695              : ! **************************************************************************************************
    1696              : !> \brief ...
    1697              : !> \param cosmat ...
    1698              : !> \param sinmat ...
    1699              : !> \param mos ...
    1700              : !> \param op_fm_set ...
    1701              : !> \param opvec ...
    1702              : ! **************************************************************************************************
    1703         1386 :    SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
    1704              : 
    1705              :       TYPE(dbcsr_type), POINTER                          :: cosmat, sinmat
    1706              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1707              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: op_fm_set
    1708              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: opvec
    1709              : 
    1710              :       INTEGER                                            :: i, nao, nmo
    1711              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1712              : 
    1713         2814 :       DO i = 1, SIZE(op_fm_set, 2) ! spin
    1714         1428 :          CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
    1715         1428 :          CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
    1716              :          CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
    1717         1428 :                             op_fm_set(1, i))
    1718         1428 :          CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
    1719              :          CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
    1720         4242 :                             op_fm_set(2, i))
    1721              :       END DO
    1722              : 
    1723         1386 :    END SUBROUTINE op_orbbas
    1724              : 
    1725              : ! **************************************************************************************************
    1726              : !> \brief ...
    1727              : !> \param cosmat ...
    1728              : !> \param sinmat ...
    1729              : !> \param mos ...
    1730              : !> \param op_fm_set ...
    1731              : !> \param mos_new ...
    1732              : ! **************************************************************************************************
    1733           48 :    SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
    1734              : 
    1735              :       TYPE(dbcsr_type), POINTER                          :: cosmat, sinmat
    1736              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1737              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: op_fm_set
    1738              :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
    1739              : 
    1740              :       INTEGER                                            :: i, icol, lcol, nao, newdim, nmo
    1741              :       LOGICAL                                            :: double_col, double_row
    1742              :       TYPE(cp_fm_struct_type), POINTER                   :: newstruct, newstruct1
    1743              :       TYPE(cp_fm_type)                                   :: work, work1, work2
    1744              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1745              : 
    1746          120 :       DO i = 1, SIZE(op_fm_set, 2) ! spin
    1747           72 :          CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
    1748           72 :          CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
    1749           72 :          double_col = .TRUE.
    1750           72 :          double_row = .FALSE.
    1751              :          CALL cp_fm_struct_double(newstruct, &
    1752              :                                   mos_new(2*i)%matrix_struct, &
    1753              :                                   mos_new(2*i)%matrix_struct%context, &
    1754              :                                   double_col, &
    1755           72 :                                   double_row)
    1756              : 
    1757           72 :          CALL cp_fm_create(work, matrix_struct=newstruct)
    1758           72 :          CALL cp_fm_create(work1, matrix_struct=newstruct)
    1759           72 :          CALL cp_fm_create(work2, matrix_struct=newstruct)
    1760           72 :          CALL cp_fm_get_info(work, ncol_global=newdim)
    1761              : 
    1762           72 :          CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
    1763          336 :          DO icol = 1, lcol
    1764         3300 :             work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
    1765         3372 :             work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
    1766              :          END DO
    1767              : 
    1768           72 :          CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
    1769           72 :          CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
    1770              : 
    1771          336 :          DO icol = 1, lcol
    1772         3300 :             work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
    1773         3372 :             work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
    1774              :          END DO
    1775              : 
    1776           72 :          CALL cp_fm_release(work1)
    1777           72 :          CALL cp_fm_release(work2)
    1778              : 
    1779              :          CALL cp_fm_struct_double(newstruct1, &
    1780              :                                   op_fm_set(1, i)%matrix_struct, &
    1781              :                                   op_fm_set(1, i)%matrix_struct%context, &
    1782              :                                   double_col, &
    1783           72 :                                   double_row)
    1784              : 
    1785           72 :          CALL cp_fm_create(work1, matrix_struct=newstruct1)
    1786              : 
    1787              :          CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
    1788           72 :                             work, 0.0_dp, work1)
    1789              : 
    1790          336 :          DO icol = 1, lcol
    1791          756 :             op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
    1792          828 :             op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
    1793              :          END DO
    1794              : 
    1795              :          CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
    1796           72 :                             work, 0.0_dp, work1)
    1797              : 
    1798          336 :          DO icol = 1, lcol
    1799              :             op_fm_set(1, i)%local_data(:, icol) = &
    1800          756 :                op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
    1801              :             op_fm_set(2, i)%local_data(:, icol) = &
    1802          828 :                op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
    1803              :          END DO
    1804              : 
    1805           72 :          CALL cp_fm_release(work)
    1806           72 :          CALL cp_fm_release(work1)
    1807           72 :          CALL cp_fm_struct_release(newstruct)
    1808          336 :          CALL cp_fm_struct_release(newstruct1)
    1809              : 
    1810              :       END DO
    1811              : 
    1812           48 :    END SUBROUTINE op_orbbas_rtp
    1813              : 
    1814              : ! **************************************************************************************************
    1815              : !> \brief ...
    1816              : !> \param qs_env ...
    1817              : !> \param magnetic ...
    1818              : !> \param nmoments ...
    1819              : !> \param reference ...
    1820              : !> \param ref_point ...
    1821              : !> \param unit_number ...
    1822              : !> \param vel_reprs ...
    1823              : !> \param com_nl ...
    1824              : ! **************************************************************************************************
    1825         1082 :    SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
    1826              : 
    1827              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1828              :       LOGICAL, INTENT(IN)                                :: magnetic
    1829              :       INTEGER, INTENT(IN)                                :: nmoments, reference
    1830              :       REAL(dp), DIMENSION(:), INTENT(IN), POINTER        :: ref_point
    1831              :       INTEGER, INTENT(IN)                                :: unit_number
    1832              :       LOGICAL, INTENT(IN), OPTIONAL                      :: vel_reprs, com_nl
    1833              : 
    1834              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_moment_locop'
    1835              : 
    1836         1082 :       CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:)        :: rlab
    1837              :       CHARACTER(LEN=default_string_length)               :: description
    1838              :       INTEGER                                            :: akind, handle, i, ia, iatom, idir, &
    1839              :                                                             ikind, ispin, ix, iy, iz, l, nm, nmom, &
    1840              :                                                             order
    1841              :       LOGICAL                                            :: my_com_nl, my_velreprs
    1842              :       REAL(dp)                                           :: charge, dd, strace, trace
    1843         1082 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
    1844         1082 :                                                             nlcom_rv, nlcom_rvr, nlcom_rxrv, &
    1845         1082 :                                                             qupole_der, rmom_vel
    1846         1082 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: rmom
    1847              :       REAL(dp), DIMENSION(3)                             :: rcc, ria
    1848              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1849              :       TYPE(cell_type), POINTER                           :: cell
    1850              :       TYPE(cp_result_type), POINTER                      :: results
    1851         1082 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: magmom, matrix_s, moments, momentum, &
    1852         1082 :                                                             rho_ao
    1853         1082 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: moments_der
    1854              :       TYPE(dbcsr_type), POINTER                          :: tmp_ao
    1855              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1856              :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1857              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1858              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1859         1082 :          POINTER                                         :: sab_all, sab_orb
    1860         1082 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1861         1082 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1862              :       TYPE(qs_rho_type), POINTER                         :: rho
    1863              : 
    1864            0 :       CPASSERT(ASSOCIATED(qs_env))
    1865              : 
    1866         1082 :       CALL timeset(routineN, handle)
    1867              : 
    1868         1082 :       my_velreprs = .FALSE.
    1869         1082 :       IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
    1870         1082 :       IF (PRESENT(com_nl)) my_com_nl = com_nl
    1871         1082 :       IF (my_velreprs) CALL cite_reference(Mattiat2019)
    1872              : 
    1873              :       ! reference point
    1874         1082 :       CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
    1875              : 
    1876              :       ! only allow for moments up to maxl set by basis
    1877         1082 :       nmom = MIN(nmoments, current_maxl)
    1878              :       ! electronic contribution
    1879         1082 :       NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
    1880              :       CALL get_qs_env(qs_env, &
    1881              :                       dft_control=dft_control, &
    1882              :                       rho=rho, &
    1883              :                       cell=cell, &
    1884              :                       results=results, &
    1885              :                       particle_set=particle_set, &
    1886              :                       qs_kind_set=qs_kind_set, &
    1887              :                       para_env=para_env, &
    1888              :                       matrix_s=matrix_s, &
    1889              :                       sab_all=sab_all, &
    1890         1082 :                       sab_orb=sab_orb)
    1891              : 
    1892         1082 :       IF (my_com_nl) THEN
    1893           40 :          IF ((nmom >= 1) .AND. my_velreprs) THEN
    1894           40 :             ALLOCATE (nlcom_rv(3))
    1895           40 :             nlcom_rv(:) = 0._dp
    1896              :          END IF
    1897           40 :          IF ((nmom >= 2) .AND. my_velreprs) THEN
    1898           40 :             ALLOCATE (nlcom_rrv(6))
    1899           40 :             nlcom_rrv(:) = 0._dp
    1900           40 :             ALLOCATE (nlcom_rvr(6))
    1901           40 :             nlcom_rvr(:) = 0._dp
    1902           40 :             ALLOCATE (nlcom_rrv_vrr(6))
    1903           40 :             nlcom_rrv_vrr(:) = 0._dp
    1904              :          END IF
    1905           40 :          IF (magnetic) THEN
    1906           18 :             ALLOCATE (nlcom_rxrv(3))
    1907           18 :             nlcom_rxrv = 0._dp
    1908              :          END IF
    1909              :          ! Calculate non local correction terms
    1910           40 :          CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
    1911              :       END IF
    1912              : 
    1913         1082 :       NULLIFY (moments)
    1914         1082 :       nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
    1915         1082 :       CALL dbcsr_allocate_matrix_set(moments, nm)
    1916         4630 :       DO i = 1, nm
    1917         3548 :          ALLOCATE (moments(i)%matrix)
    1918         3548 :          IF (my_velreprs .AND. (nmom >= 2)) THEN
    1919              :             CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
    1920          360 :                               matrix_type=dbcsr_type_symmetric)
    1921          360 :             CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
    1922              :          ELSE
    1923         3188 :             CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
    1924              :          END IF
    1925         4630 :          CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
    1926              :       END DO
    1927              : 
    1928              :       ! calculate derivatives if quadrupole in vel. reprs. is requested
    1929         1082 :       IF (my_velreprs .AND. (nmom >= 2)) THEN
    1930           40 :          NULLIFY (moments_der)
    1931           40 :          CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
    1932          160 :          DO i = 1, 3 ! x, y, z
    1933          520 :             DO idir = 1, 3 ! d/dx, d/dy, d/dz
    1934          360 :                CALL dbcsr_init_p(moments_der(i, idir)%matrix)
    1935              :                CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
    1936          360 :                                  matrix_type=dbcsr_type_antisymmetric)
    1937          360 :                CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
    1938          480 :                CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
    1939              :             END DO
    1940              :          END DO
    1941           40 :          CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
    1942              :       ELSE
    1943         1042 :          CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
    1944              :       END IF
    1945              : 
    1946         1082 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
    1947              : 
    1948         4328 :       ALLOCATE (rmom(nm + 1, 3))
    1949         3246 :       ALLOCATE (rlab(nm + 1))
    1950         1082 :       rmom = 0.0_dp
    1951         5712 :       rlab = ""
    1952              : 
    1953         1082 :       IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
    1954              :          ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
    1955              :          ! symmetric matrices)
    1956           42 :          NULLIFY (tmp_ao)
    1957           42 :          CALL dbcsr_init_p(tmp_ao)
    1958           42 :          CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
    1959           42 :          CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
    1960           42 :          CALL dbcsr_set(tmp_ao, 0.0_dp)
    1961              :       END IF
    1962              : 
    1963         1082 :       trace = 0.0_dp
    1964         2242 :       DO ispin = 1, dft_control%nspins
    1965         1160 :          CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
    1966         2242 :          rmom(1, 1) = rmom(1, 1) + trace
    1967              :       END DO
    1968              : 
    1969         4630 :       DO i = 1, SIZE(moments)
    1970         3548 :          strace = 0._dp
    1971         7330 :          DO ispin = 1, dft_control%nspins
    1972         3782 :             IF (my_velreprs .AND. nmoments >= 2) THEN
    1973              :                CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
    1974          360 :                                    0.0_dp, tmp_ao)
    1975          360 :                CALL dbcsr_trace(tmp_ao, trace)
    1976              :             ELSE
    1977         3422 :                CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
    1978              :             END IF
    1979         7330 :             strace = strace + trace
    1980              :          END DO
    1981         4630 :          rmom(i + 1, 1) = strace
    1982              :       END DO
    1983              : 
    1984         1082 :       CALL dbcsr_deallocate_matrix_set(moments)
    1985              : 
    1986              :       ! nuclear contribution
    1987              :       CALL get_qs_env(qs_env=qs_env, &
    1988         1082 :                       local_particles=local_particles)
    1989         3400 :       DO ikind = 1, SIZE(local_particles%n_el)
    1990         5178 :          DO ia = 1, local_particles%n_el(ikind)
    1991         1778 :             iatom = local_particles%list(ikind)%array(ia)
    1992              :             ! fold atomic positions back into unit cell
    1993        14224 :             ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
    1994         7112 :             ria = ria - rcc
    1995         1778 :             atomic_kind => particle_set(iatom)%atomic_kind
    1996         1778 :             CALL get_atomic_kind(atomic_kind, kind_number=akind)
    1997         1778 :             CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
    1998         1778 :             rmom(1, 2) = rmom(1, 2) - charge
    1999        10393 :             DO l = 1, nm
    2000         6297 :                ix = indco(1, l + 1)
    2001         6297 :                iy = indco(2, l + 1)
    2002         6297 :                iz = indco(3, l + 1)
    2003         6297 :                dd = 1._dp
    2004         6297 :                IF (ix > 0) dd = dd*ria(1)**ix
    2005         6297 :                IF (iy > 0) dd = dd*ria(2)**iy
    2006         6297 :                IF (iz > 0) dd = dd*ria(3)**iz
    2007         6297 :                rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
    2008         8075 :                CALL set_label(rlab(l + 1), ix, iy, iz)
    2009              :             END DO
    2010              :          END DO
    2011              :       END DO
    2012         1082 :       CALL para_env%sum(rmom(:, 2))
    2013        18218 :       rmom(:, :) = -rmom(:, :)
    2014         5712 :       rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
    2015              : 
    2016              :       ! magnetic moments
    2017         1082 :       IF (magnetic) THEN
    2018           20 :          NULLIFY (magmom)
    2019           20 :          CALL dbcsr_allocate_matrix_set(magmom, 3)
    2020           80 :          DO i = 1, SIZE(magmom)
    2021           60 :             CALL dbcsr_init_p(magmom(i)%matrix)
    2022              :             CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
    2023           60 :                               matrix_type=dbcsr_type_antisymmetric)
    2024           60 :             CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
    2025           80 :             CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
    2026              :          END DO
    2027              : 
    2028           20 :          CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
    2029              : 
    2030           60 :          ALLOCATE (mmom(SIZE(magmom)))
    2031           20 :          mmom(:) = 0.0_dp
    2032           20 :          IF (qs_env%run_rtp) THEN
    2033              :             ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
    2034              :             ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
    2035              :             ! There may be other cases, where the imaginary part of the density is relevant
    2036           12 :             NULLIFY (rho_ao)
    2037           12 :             CALL qs_rho_get(rho, rho_ao_im=rho_ao)
    2038              :          END IF
    2039              :          ! if the density is purely real this is an expensive way to calculate zero
    2040           80 :          DO i = 1, SIZE(magmom)
    2041           60 :             strace = 0._dp
    2042          120 :             DO ispin = 1, dft_control%nspins
    2043           60 :                CALL dbcsr_set(tmp_ao, 0.0_dp)
    2044              :                CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
    2045           60 :                                    0.0_dp, tmp_ao)
    2046           60 :                CALL dbcsr_trace(tmp_ao, trace)
    2047          120 :                strace = strace + trace
    2048              :             END DO
    2049           80 :             mmom(i) = strace
    2050              :          END DO
    2051              : 
    2052           20 :          CALL dbcsr_deallocate_matrix_set(magmom)
    2053              :       END IF
    2054              : 
    2055              :       ! velocity representations
    2056         1082 :       IF (my_velreprs) THEN
    2057          120 :          ALLOCATE (rmom_vel(nm))
    2058           40 :          rmom_vel = 0.0_dp
    2059              : 
    2060          120 :          DO order = 1, nmom
    2061           40 :             SELECT CASE (order)
    2062              : 
    2063              :             CASE (1) ! expectation value of momentum
    2064           40 :                NULLIFY (momentum)
    2065           40 :                CALL dbcsr_allocate_matrix_set(momentum, 3)
    2066          160 :                DO i = 1, 3
    2067          120 :                   CALL dbcsr_init_p(momentum(i)%matrix)
    2068              :                   CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
    2069          120 :                                     matrix_type=dbcsr_type_antisymmetric)
    2070          120 :                   CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
    2071          160 :                   CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
    2072              :                END DO
    2073           40 :                CALL build_lin_mom_matrix(qs_env, momentum)
    2074              : 
    2075              :                ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
    2076           40 :                IF (qs_env%run_rtp) THEN
    2077           30 :                   NULLIFY (rho_ao)
    2078           30 :                   CALL qs_rho_get(rho, rho_ao_im=rho_ao)
    2079          120 :                   DO idir = 1, SIZE(momentum)
    2080           90 :                      strace = 0._dp
    2081          180 :                      DO ispin = 1, dft_control%nspins
    2082           90 :                         CALL dbcsr_set(tmp_ao, 0.0_dp)
    2083              :                         CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
    2084           90 :                                             0.0_dp, tmp_ao)
    2085           90 :                         CALL dbcsr_trace(tmp_ao, trace)
    2086          180 :                         strace = strace + trace
    2087              :                      END DO
    2088          120 :                      rmom_vel(idir) = rmom_vel(idir) + strace
    2089              :                   END DO
    2090              :                END IF
    2091              : 
    2092           40 :                CALL dbcsr_deallocate_matrix_set(momentum)
    2093              : 
    2094              :             CASE (2) ! expectation value of quadrupole moment in vel. reprs.
    2095           40 :                ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
    2096           40 :                qupole_der = 0._dp
    2097              : 
    2098           40 :                NULLIFY (rho_ao)
    2099           40 :                CALL qs_rho_get(rho, rho_ao=rho_ao)
    2100              : 
    2101              :                ! Calculate expectation value over real part
    2102           40 :                trace = 0._dp
    2103          160 :                DO i = 1, 3
    2104          520 :                   DO idir = 1, 3
    2105          360 :                      strace = 0._dp
    2106          720 :                      DO ispin = 1, dft_control%nspins
    2107          360 :                         CALL dbcsr_set(tmp_ao, 0._dp)
    2108          360 :                         CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
    2109          360 :                         CALL dbcsr_trace(tmp_ao, trace)
    2110          720 :                         strace = strace + trace
    2111              :                      END DO
    2112          480 :                      qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
    2113              :                   END DO
    2114              :                END DO
    2115              : 
    2116           40 :                IF (qs_env%run_rtp) THEN
    2117           30 :                   NULLIFY (rho_ao)
    2118           30 :                   CALL qs_rho_get(rho, rho_ao_im=rho_ao)
    2119              : 
    2120              :                   ! Calculate expectation value over imaginary part
    2121           30 :                   trace = 0._dp
    2122          120 :                   DO i = 1, 3
    2123          390 :                      DO idir = 1, 3
    2124          270 :                         strace = 0._dp
    2125          540 :                         DO ispin = 1, dft_control%nspins
    2126          270 :                            CALL dbcsr_set(tmp_ao, 0._dp)
    2127          270 :                            CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
    2128          270 :                            CALL dbcsr_trace(tmp_ao, trace)
    2129          540 :                            strace = strace + trace
    2130              :                         END DO
    2131          360 :                         qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
    2132              :                      END DO
    2133              :                   END DO
    2134              :                END IF
    2135              : 
    2136              :                ! calculate vel. reprs. of quadrupole moment from derivatives
    2137           40 :                rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
    2138           40 :                rmom_vel(5) = -qupole_der(2) - qupole_der(4)
    2139           40 :                rmom_vel(6) = -qupole_der(3) - qupole_der(7)
    2140           40 :                rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
    2141           40 :                rmom_vel(8) = -qupole_der(6) - qupole_der(8)
    2142           40 :                rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
    2143              : 
    2144          120 :                DEALLOCATE (qupole_der)
    2145              :             CASE DEFAULT
    2146              :             END SELECT
    2147              :          END DO
    2148              :       END IF
    2149              : 
    2150         1082 :       IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
    2151           42 :          CALL dbcsr_deallocate_matrix(tmp_ao)
    2152              :       END IF
    2153         1082 :       IF (my_velreprs .AND. (nmoments >= 2)) THEN
    2154           40 :          CALL dbcsr_deallocate_matrix_set(moments_der)
    2155              :       END IF
    2156              : 
    2157         1082 :       description = "[DIPOLE]"
    2158         1082 :       CALL cp_results_erase(results=results, description=description)
    2159              :       CALL put_results(results=results, description=description, &
    2160         1082 :                        values=rmom(2:4, 3))
    2161              : 
    2162         1082 :       IF (magnetic .AND. my_velreprs) THEN
    2163           18 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom, rmom_vel=rmom_vel)
    2164         1064 :       ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
    2165            2 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., mmom=mmom)
    2166         1062 :       ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
    2167           22 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., rmom_vel=rmom_vel)
    2168              :       ELSE
    2169         1040 :          CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.FALSE.)
    2170              :       END IF
    2171              : 
    2172         1082 :       IF (my_com_nl) THEN
    2173           40 :          IF (magnetic) THEN
    2174           72 :             mmom(:) = nlcom_rxrv(:)
    2175              :          END IF
    2176           40 :          IF (my_velreprs .AND. (nmom >= 1)) THEN
    2177           40 :             DEALLOCATE (rmom_vel)
    2178           40 :             ALLOCATE (rmom_vel(21))
    2179          160 :             rmom_vel(1:3) = nlcom_rv
    2180              :          END IF
    2181           40 :          IF (my_velreprs .AND. (nmom >= 2)) THEN
    2182          280 :             rmom_vel(4:9) = nlcom_rrv
    2183          280 :             rmom_vel(10:15) = nlcom_rvr
    2184          280 :             rmom_vel(16:21) = nlcom_rrv_vrr
    2185              :          END IF
    2186           40 :          IF (magnetic .AND. .NOT. my_velreprs) THEN
    2187            0 :             CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
    2188           40 :          ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
    2189           22 :             CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
    2190           18 :          ELSE IF (my_velreprs .AND. magnetic) THEN
    2191           18 :             CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
    2192              :          END IF
    2193              : 
    2194              :       END IF
    2195              : 
    2196              :       IF (my_com_nl) THEN
    2197           40 :          IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
    2198           40 :          IF (nmom >= 2 .AND. my_velreprs) THEN
    2199           40 :             DEALLOCATE (nlcom_rrv)
    2200           40 :             DEALLOCATE (nlcom_rvr)
    2201           40 :             DEALLOCATE (nlcom_rrv_vrr)
    2202              :          END IF
    2203           40 :          IF (magnetic) DEALLOCATE (nlcom_rxrv)
    2204              :       END IF
    2205              : 
    2206         1082 :       DEALLOCATE (rmom)
    2207         1082 :       DEALLOCATE (rlab)
    2208         1082 :       IF (magnetic) THEN
    2209           20 :          DEALLOCATE (mmom)
    2210              :       END IF
    2211         1082 :       IF (my_velreprs) THEN
    2212           40 :          DEALLOCATE (rmom_vel)
    2213              :       END IF
    2214              : 
    2215         1082 :       CALL timestop(handle)
    2216              : 
    2217         2164 :    END SUBROUTINE qs_moment_locop
    2218              : 
    2219              : ! **************************************************************************************************
    2220              : !> \brief ...
    2221              : !> \param label ...
    2222              : !> \param ix ...
    2223              : !> \param iy ...
    2224              : !> \param iz ...
    2225              : ! **************************************************************************************************
    2226         8415 :    SUBROUTINE set_label(label, ix, iy, iz)
    2227              :       CHARACTER(LEN=*), INTENT(OUT)                      :: label
    2228              :       INTEGER, INTENT(IN)                                :: ix, iy, iz
    2229              : 
    2230              :       INTEGER                                            :: i
    2231              : 
    2232         8415 :       label = ""
    2233        11581 :       DO i = 1, ix
    2234        11581 :          WRITE (label(i:), "(A1)") "X"
    2235              :       END DO
    2236        11581 :       DO i = ix + 1, ix + iy
    2237        11581 :          WRITE (label(i:), "(A1)") "Y"
    2238              :       END DO
    2239        11581 :       DO i = ix + iy + 1, ix + iy + iz
    2240        11581 :          WRITE (label(i:), "(A1)") "Z"
    2241              :       END DO
    2242              : 
    2243         8415 :    END SUBROUTINE set_label
    2244              : 
    2245              : ! **************************************************************************************************
    2246              : !> \brief ...
    2247              : !> \param unit_number ...
    2248              : !> \param nmom ...
    2249              : !> \param rmom ...
    2250              : !> \param rlab ...
    2251              : !> \param rcc ...
    2252              : !> \param cell ...
    2253              : !> \param periodic ...
    2254              : !> \param mmom ...
    2255              : !> \param rmom_vel ...
    2256              : ! **************************************************************************************************
    2257         1578 :    SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
    2258              :       INTEGER, INTENT(IN)                                :: unit_number, nmom
    2259              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: rmom
    2260              :       CHARACTER(LEN=8), DIMENSION(:)                     :: rlab
    2261              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rcc
    2262              :       TYPE(cell_type), POINTER                           :: cell
    2263              :       LOGICAL                                            :: periodic
    2264              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: mmom, rmom_vel
    2265              : 
    2266              :       INTEGER                                            :: i, i0, i1, j, l
    2267              :       REAL(dp)                                           :: dd
    2268              : 
    2269         1578 :       IF (unit_number > 0) THEN
    2270         2438 :          DO l = 0, nmom
    2271          805 :             SELECT CASE (l)
    2272              :             CASE (0)
    2273          805 :                WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
    2274          805 :                WRITE (unit_number, "(T3,A)") "Charges"
    2275              :                WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
    2276          805 :                   "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
    2277              :             CASE (1)
    2278          805 :                IF (periodic) THEN
    2279          246 :                   WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
    2280          246 :                   WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
    2281          984 :                   WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
    2282          984 :                   WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
    2283          984 :                   WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
    2284              :                ELSE
    2285          559 :                   WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
    2286              :                END IF
    2287         3220 :                dd = SQRT(SUM(rmom(2:4, 3)**2))*debye
    2288          805 :                WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
    2289              :                WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
    2290         3220 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
    2291              :             CASE (2)
    2292           21 :                WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
    2293              :                WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2294           84 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
    2295              :                WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2296           84 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
    2297              :             CASE (3)
    2298            1 :                WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
    2299              :                WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
    2300            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
    2301              :                WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
    2302            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
    2303              :                WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
    2304            3 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
    2305              :             CASE (4)
    2306            1 :                WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
    2307              :                WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
    2308            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
    2309              :                WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
    2310            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
    2311              :                WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
    2312            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
    2313              :                WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
    2314            5 :                   (TRIM(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
    2315              :             CASE DEFAULT
    2316            0 :                WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
    2317            0 :                   "  L=", l
    2318            0 :                i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
    2319            0 :                i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
    2320            0 :                dd = debye/(bohr)**(l - 1)
    2321         1633 :                DO i = i0, i1, 3
    2322              :                   WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
    2323            0 :                      (TRIM(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, MIN(i1, i + 2))
    2324              :                END DO
    2325              :             END SELECT
    2326              :          END DO
    2327          805 :          IF (PRESENT(mmom)) THEN
    2328           28 :             IF (nmom >= 1) THEN
    2329          112 :                dd = SQRT(SUM(mmom(1:3)**2))
    2330           28 :                WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
    2331              :                WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
    2332          112 :                   (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
    2333              :             END IF
    2334              :          END IF
    2335          805 :          IF (PRESENT(rmom_vel)) THEN
    2336           96 :             DO l = 1, nmom
    2337           38 :                SELECT CASE (l)
    2338              :                CASE (1)
    2339          152 :                   dd = SQRT(SUM(rmom_vel(1:3)**2))
    2340           38 :                   WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
    2341              :                   WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
    2342          152 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
    2343              :                CASE (2)
    2344           20 :                   WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
    2345              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2346           80 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
    2347              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2348          138 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
    2349              :                CASE DEFAULT
    2350              :                END SELECT
    2351              :             END DO
    2352              :          END IF
    2353              :       END IF
    2354              : 
    2355         1578 :    END SUBROUTINE print_moments
    2356              : 
    2357              : ! **************************************************************************************************
    2358              : !> \brief ...
    2359              : !> \param unit_number ...
    2360              : !> \param nmom ...
    2361              : !> \param rlab ...
    2362              : !> \param mmom ...
    2363              : !> \param rmom_vel ...
    2364              : ! **************************************************************************************************
    2365           58 :    SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
    2366              :       INTEGER, INTENT(IN)                                :: unit_number, nmom
    2367              :       CHARACTER(LEN=8), DIMENSION(:)                     :: rlab
    2368              :       REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL       :: mmom, rmom_vel
    2369              : 
    2370              :       INTEGER                                            :: i, l
    2371              :       REAL(dp)                                           :: dd
    2372              : 
    2373           58 :       IF (unit_number > 0) THEN
    2374           38 :          IF (PRESENT(mmom)) THEN
    2375           27 :             IF (nmom >= 1) THEN
    2376          108 :                dd = SQRT(SUM(mmom(1:3)**2))
    2377           27 :                WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
    2378              :                WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
    2379          108 :                   (TRIM(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
    2380              :             END IF
    2381              :          END IF
    2382           38 :          IF (PRESENT(rmom_vel)) THEN
    2383           96 :             DO l = 1, nmom
    2384           38 :                SELECT CASE (l)
    2385              :                CASE (1)
    2386          152 :                   dd = SQRT(SUM(rmom_vel(1:3)**2))
    2387           38 :                   WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
    2388              :                   WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
    2389          152 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
    2390              :                CASE (2)
    2391           20 :                   WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
    2392              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2393           80 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
    2394              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2395           80 :                      (TRIM(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
    2396           20 :                   WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
    2397              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2398           80 :                      (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
    2399              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2400           80 :                      (TRIM(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
    2401           20 :                   WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
    2402              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2403           80 :                      (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
    2404              :                   WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
    2405          138 :                      (TRIM(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
    2406              :                CASE DEFAULT
    2407              :                END SELECT
    2408              :             END DO
    2409              :          END IF
    2410              :       END IF
    2411              : 
    2412           58 :    END SUBROUTINE print_moments_nl
    2413              : 
    2414              : ! **************************************************************************************************
    2415              : !> \brief Calculate the expectation value of operators related to non-local potential:
    2416              : !>        [r, Vnl], noted rv
    2417              : !>        r x [r,Vnl], noted rxrv
    2418              : !>        [rr,Vnl], noted rrv
    2419              : !>        r x Vnl x r, noted rvr
    2420              : !>        r x r x Vnl + Vnl x r x r, noted rrv_vrr
    2421              : !>        Note that the 3 first operator are commutator while the 2 last
    2422              : !>        are not. For reading clarity the same notation is used for all 5
    2423              : !>        operators.
    2424              : !> \param qs_env ...
    2425              : !> \param nlcom_rv ...
    2426              : !> \param nlcom_rxrv ...
    2427              : !> \param nlcom_rrv ...
    2428              : !> \param nlcom_rvr ...
    2429              : !> \param nlcom_rrv_vrr ...
    2430              : !> \param ref_point ...
    2431              : ! **************************************************************************************************
    2432           76 :    SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
    2433              :                                             nlcom_rrv_vrr, ref_point)
    2434              : 
    2435              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2436              :       REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL      :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
    2437              :                                                             nlcom_rvr, nlcom_rrv_vrr
    2438              :       REAL(dp), DIMENSION(3)                             :: ref_point
    2439              : 
    2440              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_commutator_nl_terms'
    2441              : 
    2442              :       INTEGER                                            :: handle, ind, ispin
    2443              :       LOGICAL                                            :: calc_rrv, calc_rrv_vrr, calc_rv, &
    2444              :                                                             calc_rvr, calc_rxrv
    2445              :       REAL(dp)                                           :: eps_ppnl, strace, trace
    2446              :       TYPE(cell_type), POINTER                           :: cell
    2447           76 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
    2448           76 :                                                             matrix_rvr, matrix_rxrv, matrix_s, &
    2449           76 :                                                             rho_ao
    2450              :       TYPE(dbcsr_type), POINTER                          :: tmp_ao
    2451              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2452              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2453           76 :          POINTER                                         :: sab_all, sab_orb, sap_ppnl
    2454           76 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2455           76 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2456              :       TYPE(qs_rho_type), POINTER                         :: rho
    2457              : 
    2458           76 :       CALL timeset(routineN, handle)
    2459              : 
    2460           76 :       calc_rv = .FALSE.
    2461           76 :       calc_rxrv = .FALSE.
    2462           76 :       calc_rrv = .FALSE.
    2463           76 :       calc_rvr = .FALSE.
    2464           76 :       calc_rrv_vrr = .FALSE.
    2465              : 
    2466              :       ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
    2467              :       ! The real part of the density matrix rho_ao is symmetric so that
    2468              :       ! the expectation value of real density matrix is zero. Hence, if
    2469              :       ! the density matrix is real, no need to compute these quantities.
    2470              :       ! This is not the case for rvr and rrv_vrr which are symmetric.
    2471              : 
    2472           76 :       IF (ALLOCATED(nlcom_rv)) THEN
    2473           76 :          nlcom_rv(:) = 0._dp
    2474           76 :          IF (qs_env%run_rtp) calc_rv = .TRUE.
    2475              :       END IF
    2476           76 :       IF (ALLOCATED(nlcom_rxrv)) THEN
    2477           54 :          nlcom_rxrv(:) = 0._dp
    2478           54 :          IF (qs_env%run_rtp) calc_rxrv = .TRUE.
    2479              :       END IF
    2480           76 :       IF (ALLOCATED(nlcom_rrv)) THEN
    2481           40 :          nlcom_rrv(:) = 0._dp
    2482           40 :          IF (qs_env%run_rtp) calc_rrv = .TRUE.
    2483              :       END IF
    2484           76 :       IF (ALLOCATED(nlcom_rvr)) THEN
    2485           40 :          nlcom_rvr(:) = 0._dp
    2486           40 :          calc_rvr = .TRUE.
    2487              :       END IF
    2488           76 :       IF (ALLOCATED(nlcom_rrv_vrr)) THEN
    2489           40 :          nlcom_rrv_vrr(:) = 0._dp
    2490           40 :          calc_rrv_vrr = .TRUE.
    2491              :       END IF
    2492              : 
    2493           76 :       IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv .OR. calc_rvr .OR. calc_rrv_vrr)) THEN
    2494           12 :          CALL timestop(handle)
    2495           12 :          RETURN
    2496              :       END IF
    2497              : 
    2498           64 :       NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
    2499              :       CALL get_qs_env(qs_env, &
    2500              :                       cell=cell, &
    2501              :                       dft_control=dft_control, &
    2502              :                       matrix_s=matrix_s, &
    2503              :                       particle_set=particle_set, &
    2504              :                       qs_kind_set=qs_kind_set, &
    2505              :                       rho=rho, &
    2506              :                       sab_orb=sab_orb, &
    2507              :                       sab_all=sab_all, &
    2508           64 :                       sap_ppnl=sap_ppnl)
    2509              : 
    2510           64 :       eps_ppnl = dft_control%qs_control%eps_ppnl
    2511              : 
    2512              :       ! Allocate storage
    2513           64 :       NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
    2514           64 :       IF (calc_rv) THEN
    2515           54 :          CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
    2516          216 :          DO ind = 1, 3
    2517          162 :             CALL dbcsr_init_p(matrix_rv(ind)%matrix)
    2518              :             CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
    2519          162 :                               matrix_type=dbcsr_type_antisymmetric)
    2520          162 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
    2521          216 :             CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
    2522              :          END DO
    2523              :       END IF
    2524              : 
    2525           64 :       IF (calc_rxrv) THEN
    2526           36 :          CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
    2527          144 :          DO ind = 1, 3
    2528          108 :             CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
    2529              :             CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
    2530          108 :                               matrix_type=dbcsr_type_antisymmetric)
    2531          108 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
    2532          144 :             CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
    2533              :          END DO
    2534              :       END IF
    2535              : 
    2536           64 :       IF (calc_rrv) THEN
    2537           30 :          CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
    2538          210 :          DO ind = 1, 6
    2539          180 :             CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
    2540              :             CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
    2541          180 :                               matrix_type=dbcsr_type_antisymmetric)
    2542          180 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
    2543          210 :             CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
    2544              :          END DO
    2545              :       END IF
    2546              : 
    2547           64 :       IF (calc_rvr) THEN
    2548           40 :          CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
    2549          280 :          DO ind = 1, 6
    2550          240 :             CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
    2551              :             CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
    2552          240 :                               matrix_type=dbcsr_type_symmetric)
    2553          240 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
    2554          280 :             CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
    2555              :          END DO
    2556              :       END IF
    2557           64 :       IF (calc_rrv_vrr) THEN
    2558           40 :          CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
    2559          280 :          DO ind = 1, 6
    2560          240 :             CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
    2561              :             CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
    2562          240 :                               matrix_type=dbcsr_type_symmetric)
    2563          240 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
    2564          280 :             CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
    2565              :          END DO
    2566              :       END IF
    2567              : 
    2568              :       ! calculate evaluation of operators in AO basis set
    2569              :       CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
    2570              :                             matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
    2571           64 :                             matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
    2572              : 
    2573              :       ! Calculate expectation values
    2574              :       ! Real part
    2575           64 :       NULLIFY (tmp_ao)
    2576           64 :       CALL dbcsr_init_p(tmp_ao)
    2577           64 :       CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
    2578           64 :       CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
    2579           64 :       CALL dbcsr_set(tmp_ao, 0.0_dp)
    2580              : 
    2581           64 :       IF (calc_rvr .OR. calc_rrv_vrr) THEN
    2582           40 :          NULLIFY (rho_ao)
    2583           40 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
    2584              : 
    2585           40 :          IF (calc_rvr) THEN
    2586              :             trace = 0._dp
    2587          280 :             DO ind = 1, SIZE(matrix_rvr)
    2588          240 :                strace = 0._dp
    2589          480 :                DO ispin = 1, dft_control%nspins
    2590          240 :                   CALL dbcsr_set(tmp_ao, 0.0_dp)
    2591              :                   CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
    2592          240 :                                       0.0_dp, tmp_ao)
    2593          240 :                   CALL dbcsr_trace(tmp_ao, trace)
    2594          480 :                   strace = strace + trace
    2595              :                END DO
    2596          280 :                nlcom_rvr(ind) = nlcom_rvr(ind) + strace
    2597              :             END DO
    2598              :          END IF
    2599              : 
    2600           40 :          IF (calc_rrv_vrr) THEN
    2601              :             trace = 0._dp
    2602          280 :             DO ind = 1, SIZE(matrix_rrv_vrr)
    2603          240 :                strace = 0._dp
    2604          480 :                DO ispin = 1, dft_control%nspins
    2605          240 :                   CALL dbcsr_set(tmp_ao, 0.0_dp)
    2606              :                   CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
    2607          240 :                                       0.0_dp, tmp_ao)
    2608          240 :                   CALL dbcsr_trace(tmp_ao, trace)
    2609          480 :                   strace = strace + trace
    2610              :                END DO
    2611          280 :                nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
    2612              :             END DO
    2613              :          END IF
    2614              :       END IF
    2615              : 
    2616              :       ! imagninary part of the density matrix
    2617           64 :       NULLIFY (rho_ao)
    2618           64 :       CALL qs_rho_get(rho, rho_ao_im=rho_ao)
    2619              : 
    2620           64 :       IF (calc_rv) THEN
    2621              :          trace = 0._dp
    2622          216 :          DO ind = 1, SIZE(matrix_rv)
    2623          162 :             strace = 0._dp
    2624          324 :             DO ispin = 1, dft_control%nspins
    2625          162 :                CALL dbcsr_set(tmp_ao, 0.0_dp)
    2626              :                CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
    2627          162 :                                    0.0_dp, tmp_ao)
    2628          162 :                CALL dbcsr_trace(tmp_ao, trace)
    2629          324 :                strace = strace + trace
    2630              :             END DO
    2631          216 :             nlcom_rv(ind) = nlcom_rv(ind) + strace
    2632              :          END DO
    2633              :       END IF
    2634              : 
    2635           64 :       IF (calc_rrv) THEN
    2636              :          trace = 0._dp
    2637          210 :          DO ind = 1, SIZE(matrix_rrv)
    2638          180 :             strace = 0._dp
    2639          360 :             DO ispin = 1, dft_control%nspins
    2640          180 :                CALL dbcsr_set(tmp_ao, 0.0_dp)
    2641              :                CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
    2642          180 :                                    0.0_dp, tmp_ao)
    2643          180 :                CALL dbcsr_trace(tmp_ao, trace)
    2644          360 :                strace = strace + trace
    2645              :             END DO
    2646          210 :             nlcom_rrv(ind) = nlcom_rrv(ind) + strace
    2647              :          END DO
    2648              :       END IF
    2649              : 
    2650           64 :       IF (calc_rxrv) THEN
    2651              :          trace = 0._dp
    2652          144 :          DO ind = 1, SIZE(matrix_rxrv)
    2653          108 :             strace = 0._dp
    2654          216 :             DO ispin = 1, dft_control%nspins
    2655          108 :                CALL dbcsr_set(tmp_ao, 0.0_dp)
    2656              :                CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
    2657          108 :                                    0.0_dp, tmp_ao)
    2658          108 :                CALL dbcsr_trace(tmp_ao, trace)
    2659          216 :                strace = strace + trace
    2660              :             END DO
    2661          144 :             nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
    2662              :          END DO
    2663              :       END IF
    2664           64 :       CALL dbcsr_deallocate_matrix(tmp_ao)
    2665           64 :       IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
    2666           64 :       IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
    2667           64 :       IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
    2668           64 :       IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
    2669           64 :       IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
    2670              : 
    2671           64 :       CALL timestop(handle)
    2672           76 :    END SUBROUTINE calculate_commutator_nl_terms
    2673              : 
    2674              : ! **************************************************************************************************
    2675              : !> \brief Get list of kpoints from input to compute dipole moment elements
    2676              : !> \param qs_env ...
    2677              : !> \param xkp ...
    2678              : !> \param special_pnts ...
    2679              : !> \author Shridhar Shanbhag
    2680              : ! **************************************************************************************************
    2681           10 :    SUBROUTINE get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
    2682              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2683              :       TYPE(section_vals_type), POINTER                   :: kpnts, kpset
    2684              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: cart_hmat, hmat
    2685              : 
    2686              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_xkp_for_dipole_calc'
    2687              : 
    2688              :       CHARACTER(LEN=default_string_length)               :: ustr
    2689              :       TYPE(kpoint_type), POINTER                         :: kpoint_work
    2690              :       TYPE(cell_type), POINTER                           :: cell
    2691              :       CHARACTER(LEN=default_string_length), &
    2692           10 :          DIMENSION(:), POINTER                           :: strptr
    2693              :       CHARACTER(LEN=default_string_length), &
    2694           10 :          DIMENSION(:), POINTER                           :: special_pnts, spname
    2695              :       CHARACTER(LEN=max_line_length)                     :: error_message
    2696              :       INTEGER                                            :: handle, i, ik, ikk, ip, &
    2697              :                                                             n_ptr, npline, nkp
    2698              :       LOGICAL                                            :: explicit_kpnts, explicit_kpset
    2699           10 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: kspecial, xkp
    2700              :       REAL(KIND=dp), DIMENSION(3)                        :: kpptr
    2701           10 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2702              : 
    2703           10 :       CALL timeset(routineN, handle)
    2704           10 :       kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
    2705           10 :       kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
    2706           10 :       CALL section_vals_get(kpset, explicit=explicit_kpset)
    2707           10 :       CALL section_vals_get(kpnts, explicit=explicit_kpnts)
    2708           10 :       IF (explicit_kpset .AND. explicit_kpnts) then
    2709            0 :          CPABORT("Both KPOINT_SET and KPOINTS present in MOMENTS section")
    2710              :       end if
    2711              : 
    2712           10 :       IF (explicit_kpset) THEN
    2713            4 :          CALL get_qs_env(qs_env, cell=cell)
    2714            4 :          CALL get_cell(cell, h=hmat)
    2715            4 :          cart_hmat(:, :) = hmat(:, :)
    2716            4 :          IF (cell%input_cell_canonicalized) cart_hmat(:, :) = cell%input_hmat(:, :)
    2717            4 :          CALL section_vals_val_get(kpset, "NPOINTS", i_val=npline)
    2718            4 :          CALL section_vals_val_get(kpset, "UNITS", c_val=ustr)
    2719            4 :          CALL uppercase(ustr)
    2720            4 :          CALL section_vals_val_get(kpset, "SPECIAL_POINT", n_rep_val=n_ptr)
    2721            4 :          CPASSERT(n_ptr > 0)
    2722           12 :          ALLOCATE (kspecial(3, n_ptr))
    2723           12 :          ALLOCATE (spname(n_ptr))
    2724            8 :          DO ip = 1, n_ptr
    2725            4 :             CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_val=ip, c_vals=strptr)
    2726            4 :             IF (SIZE(strptr(:), 1) == 4) THEN
    2727            2 :                spname(ip) = strptr(1)
    2728            8 :                DO i = 1, 3
    2729            6 :                   CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
    2730            8 :                   IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
    2731              :                END DO
    2732            2 :             ELSE IF (SIZE(strptr(:), 1) == 3) THEN
    2733            2 :                spname(ip) = "not specified"
    2734            8 :                DO i = 1, 3
    2735            6 :                   CALL read_float_object(strptr(i), kpptr(i), error_message)
    2736            8 :                   IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message))
    2737              :                END DO
    2738              :             ELSE
    2739            0 :                CPABORT("Input SPECIAL_POINT invalid")
    2740              :             END IF
    2741            4 :             SELECT CASE (ustr)
    2742              :             CASE ("B_VECTOR")
    2743           16 :                kspecial(1:3, ip) = kpptr(1:3)
    2744              :             CASE ("CART_ANGSTROM")
    2745              :                kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
    2746              :                                     kpptr(2)*cart_hmat(2, 1:3) + &
    2747            0 :                                     kpptr(3)*cart_hmat(3, 1:3))/twopi*angstrom
    2748              :             CASE ("CART_BOHR")
    2749              :                kspecial(1:3, ip) = (kpptr(1)*cart_hmat(1, 1:3) + &
    2750              :                                     kpptr(2)*cart_hmat(2, 1:3) + &
    2751            0 :                                     kpptr(3)*cart_hmat(3, 1:3))/twopi
    2752              :             CASE DEFAULT
    2753            4 :                CPABORT("Unknown unit <"//TRIM(ustr)//"> specified for k-point definition")
    2754              :             END SELECT
    2755              :          END DO
    2756            4 :          nkp = (n_ptr - 1)*npline + 1
    2757            4 :          CPASSERT(nkp >= 1)
    2758              : 
    2759              :          ! Initialize environment and calculate MOs
    2760           12 :          ALLOCATE (xkp(3, nkp))
    2761           12 :          ALLOCATE (special_pnts(nkp))
    2762            8 :          special_pnts(:) = ""
    2763           16 :          xkp(1:3, 1) = kspecial(1:3, 1)
    2764            4 :          ikk = 1
    2765            4 :          special_pnts(ikk) = spname(1)
    2766            4 :          DO ik = 2, n_ptr
    2767            0 :             DO ip = 1, npline
    2768            0 :                ikk = ikk + 1
    2769              :                xkp(1:3, ikk) = kspecial(1:3, ik - 1) + &
    2770              :                                REAL(ip, KIND=dp)/REAL(npline, KIND=dp)* &
    2771            0 :                                (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
    2772              :             END DO
    2773            4 :             special_pnts(ikk) = spname(ik)
    2774              :          END DO
    2775           12 :          DEALLOCATE (spname, kspecial)
    2776            6 :       ELSE IF (explicit_kpnts) THEN
    2777            2 :          CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
    2778            2 :          CALL get_cell(cell, h=hmat)
    2779            2 :          NULLIFY (kpoint_work)
    2780            2 :          CALL kpoint_create(kpoint_work)
    2781            2 :          CALL read_kpoint_section(kpoint_work, kpnts, hmat, cell)
    2782            2 :          CALL kpoint_initialize(kpoint_work, particle_set, cell)
    2783            2 :          nkp = kpoint_work%nkp
    2784            6 :          ALLOCATE (xkp(3, nkp))
    2785            6 :          ALLOCATE (special_pnts(nkp))
    2786            4 :          special_pnts(:) = ""
    2787           10 :          xkp(1:3, :) = kpoint_work%xkp(1:3, :)
    2788            2 :          CALL kpoint_release(kpoint_work)
    2789              :       ELSE
    2790              :          ! use k-point mesh from DFT calculation
    2791            4 :          CALL get_qs_env(qs_env, kpoints=kpoint_work)
    2792            4 :          nkp = kpoint_work%nkp
    2793            4 :          nkp = kpoint_work%nkp
    2794           12 :          ALLOCATE (xkp(3, nkp))
    2795           12 :          ALLOCATE (special_pnts(nkp))
    2796          252 :          special_pnts(:) = ""
    2797          996 :          xkp(1:3, :) = kpoint_work%xkp(1:3, :)
    2798              :       END IF
    2799           10 :       CALL timestop(handle)
    2800              : 
    2801           10 :    END SUBROUTINE get_xkp_for_dipole_calc
    2802              : ! **************************************************************************************************
    2803              : !> \brief Calculate local moment matrix for a periodic system for all image cells
    2804              : !> \param qs_env ...
    2805              : !> \param moments_rs_img ...
    2806              : !> \param rcc ...
    2807              : !> \author Shridhar Shanbhag
    2808              : ! **************************************************************************************************
    2809            8 :    SUBROUTINE build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
    2810              : 
    2811              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2812              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: moments_rs_img
    2813              :       REAL(KIND=dp), DIMENSION(3), OPTIONAL              :: rcc
    2814              : 
    2815              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moment_matrix_rs_img'
    2816              : 
    2817              :       INTEGER :: handle, i_dir, iatom, ic, ikind, iset, j, jatom, jkind, jset, &
    2818              :                  ldsa, ldsb, ldwork, ncoa, ncob, nimg, nkind, nseta, nsetb, nsize, sgfa, sgfb
    2819              :       INTEGER, DIMENSION(3)                              :: icell
    2820            8 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
    2821            8 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    2822              :       LOGICAL                                            :: found
    2823              :       REAL(dp), DIMENSION(3)                             :: ra, rab, rac, rb, rbc, rc
    2824            8 :       REAL(dp), DIMENSION(:, :), POINTER                 :: dblock, work
    2825            8 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: dipab
    2826              :       REAL(KIND=dp)                                      :: dab
    2827              :       TYPE(cell_type), POINTER                           :: cell
    2828            8 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp
    2829              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2830            8 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
    2831              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set, basis_set_a, basis_set_b
    2832              :       TYPE(kpoint_type), POINTER                         :: kpoints_all
    2833              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2834              :       TYPE(neighbor_list_iterator_p_type), &
    2835            8 :          DIMENSION(:), POINTER                           :: nl_iterator
    2836              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2837            8 :          POINTER                                         :: sab_all
    2838            8 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2839            8 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2840              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    2841              : 
    2842            8 :       CALL timeset(routineN, handle)
    2843              : 
    2844              :       CALL get_qs_env(qs_env=qs_env, &
    2845              :                       dft_control=dft_control, &
    2846              :                       qs_kind_set=qs_kind_set, &
    2847              :                       matrix_ks_kp=matrix_ks_kp, &
    2848              :                       particle_set=particle_set, &
    2849              :                       cell=cell, &
    2850              :                       para_env=para_env, &
    2851            8 :                       sab_all=sab_all)
    2852              : 
    2853            8 :       NULLIFY (kpoints_all)
    2854            8 :       CALL kpoint_create(kpoints_all)
    2855            8 :       CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
    2856              : 
    2857            8 :       nkind = SIZE(qs_kind_set)
    2858           36 :       ALLOCATE (basis_set_list(nkind))
    2859           20 :       DO ikind = 1, nkind
    2860           12 :          qs_kind => qs_kind_set(ikind)
    2861           12 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
    2862           20 :          IF (ASSOCIATED(basis_set)) THEN
    2863           12 :             basis_set_list(ikind)%gto_basis_set => basis_set
    2864              :          ELSE
    2865            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
    2866              :          END IF
    2867              :       END DO
    2868              : 
    2869            8 :       rc(:) = 0._dp
    2870            8 :       IF (PRESENT(rcc)) rc(:) = rcc(:)
    2871              : 
    2872            8 :       CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list)
    2873            8 :       CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
    2874            8 :       nsize = SIZE(index_to_cell, 2)
    2875            8 :       CPASSERT(SIZE(moments_rs_img, 2) == nsize)
    2876           32 :       DO i_dir = 1, 3
    2877         2360 :          DO j = 1, nsize
    2878         2328 :             ALLOCATE (moments_rs_img(i_dir, j)%matrix)
    2879              :             CALL dbcsr_create(matrix=moments_rs_img(i_dir, j)%matrix, &
    2880              :                               template=matrix_ks_kp(1, 1)%matrix, &
    2881              :                               matrix_type=dbcsr_type_no_symmetry, &
    2882         2328 :                               name="DIPMAT")
    2883         2328 :             CALL cp_dbcsr_alloc_block_from_nbl(moments_rs_img(i_dir, j)%matrix, sab_all)
    2884         2352 :             CALL dbcsr_set(moments_rs_img(i_dir, j)%matrix, 0.0_dp)
    2885              :          END DO
    2886              :       END DO
    2887              : 
    2888            8 :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
    2889           40 :       ALLOCATE (dipab(ldwork, ldwork, 3))
    2890           32 :       ALLOCATE (work(ldwork, ldwork))
    2891              : 
    2892            8 :       CALL neighbor_list_iterator_create(nl_iterator, sab_all)
    2893         2100 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    2894              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
    2895         2092 :                                 iatom=iatom, jatom=jatom, r=rab, cell=icell)
    2896              : 
    2897         2092 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
    2898         2092 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    2899         2092 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
    2900         2092 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    2901              :          ASSOCIATE ( &
    2902              :             ! basis ikind
    2903              :             first_sgfa => basis_set_a%first_sgf, &
    2904              :             la_max => basis_set_a%lmax, &
    2905              :             la_min => basis_set_a%lmin, &
    2906              :             npgfa => basis_set_a%npgf, &
    2907              :             nsgfa => basis_set_a%nsgf_set, &
    2908              :             rpgfa => basis_set_a%pgf_radius, &
    2909              :             set_radius_a => basis_set_a%set_radius, &
    2910              :             sphi_a => basis_set_a%sphi, &
    2911              :             zeta => basis_set_a%zet, &
    2912              :             ! basis jkind, &
    2913              :             first_sgfb => basis_set_b%first_sgf, &
    2914              :             lb_max => basis_set_b%lmax, &
    2915              :             lb_min => basis_set_b%lmin, &
    2916              :             npgfb => basis_set_b%npgf, &
    2917              :             nsgfb => basis_set_b%nsgf_set, &
    2918              :             rpgfb => basis_set_b%pgf_radius, &
    2919              :             set_radius_b => basis_set_b%set_radius, &
    2920              :             sphi_b => basis_set_b%sphi, &
    2921              :             zetb => basis_set_b%zet)
    2922              : 
    2923         2092 :             nseta = basis_set_a%nset
    2924         2092 :             nsetb = basis_set_b%nset
    2925              : 
    2926         2092 :             ldsa = SIZE(sphi_a, 1)
    2927         2092 :             ldsb = SIZE(sphi_b, 1)
    2928              : 
    2929         2092 :             NULLIFY (dblock)
    2930              : 
    2931         2092 :             ra = pbc(particle_set(iatom)%r(:), cell)
    2932         8368 :             rb(:) = ra(:) + rab(:)
    2933         8368 :             rac = ra - rc
    2934         8368 :             rbc = rb - rc
    2935         8368 :             dab = norm2(rab)
    2936              : 
    2937         2092 :             ic = cell_to_index(icell(1), icell(2), icell(3))
    2938              : 
    2939         6276 :             DO iset = 1, nseta
    2940              : 
    2941         2092 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
    2942         2092 :                sgfa = first_sgfa(1, iset)
    2943              : 
    2944         6276 :                DO jset = 1, nsetb
    2945              : 
    2946         2092 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
    2947              : 
    2948         2092 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
    2949         2092 :                   sgfb = first_sgfb(1, jset)
    2950              : 
    2951              :                   CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
    2952              :                               lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), 1, &
    2953         2092 :                               rac, rbc, dipab)
    2954        10460 :                   DO i_dir = 1, 3
    2955              :                      CALL dbcsr_get_block_p(matrix=moments_rs_img(i_dir, ic)%matrix, &
    2956         6276 :                                             row=iatom, col=jatom, BLOCK=dblock, found=found)
    2957         6276 :                      CPASSERT(found)
    2958              :                      CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
    2959              :                                 1.0_dp, dipab(1, 1, i_dir), ldwork, &
    2960         6276 :                                 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
    2961              : 
    2962              :                      CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
    2963              :                                 1.0_dp, sphi_a(1, sgfa), ldsa, &
    2964        14644 :                                 work(1, 1), ldwork, 1.0_dp, dblock(1, 1), SIZE(dblock, 1))
    2965              :                   END DO
    2966              :                END DO
    2967              :             END DO
    2968              :          END ASSOCIATE
    2969              :       END DO
    2970            8 :       CALL neighbor_list_iterator_release(nl_iterator)
    2971            8 :       CALL kpoint_release(kpoints_all)
    2972            8 :       DEALLOCATE (dipab, work, basis_set_list)
    2973            8 :       CALL timestop(handle)
    2974              : 
    2975           24 :    END SUBROUTINE build_local_moment_matrix_rs_img
    2976              : 
    2977              : ! **************************************************************************************************
    2978              : !> \brief Calculates the dipole moments and berry curvature for periodic systems for kpoints
    2979              : !> \param qs_env ...
    2980              : !> \param xkp list of kpoints
    2981              : !> \param dipole ...
    2982              : !> \param rcc coordinates about which to calculate the dipole
    2983              : !> \param berry_c berry curvature calculated using Ω^γ_n = Σ_m 2*Im[d^α_nm (d^β_mn)*]
    2984              : !> \param do_parallel option to distribute the result in dipole across
    2985              : !>        different MPI ranks
    2986              : !> \author Shridhar Shanbhag
    2987              : ! **************************************************************************************************
    2988            8 :    SUBROUTINE qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
    2989              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2990              :       LOGICAL, OPTIONAL                                  :: do_parallel
    2991              :       LOGICAL                                            :: my_do_parallel, calc_bc
    2992              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_moment_kpoints_deep'
    2993              :       COMPLEX(KIND=dp)                                   :: phase, tmp_max
    2994            8 :       COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE     :: C_k, H_k, S_k, D_k, CDC, C_dH_C, &
    2995            8 :                                                             C_dS_C, dH_dk_i, dS_dk_i
    2996            8 :       COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE  :: dip
    2997              :       COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
    2998              :          ALLOCATABLE                                     :: dipole
    2999              :       INTEGER                                            :: handle, i_dir, ikp, nkp, &
    3000              :                                                             n_img_scf, n_img_all, nao, &
    3001              :                                                             num_pe, num_copy, mepos, n, m, mu, &
    3002              :                                                             ispin, nspin
    3003              :       INTEGER, DIMENSION(3)                              :: periodic
    3004            8 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell_all
    3005            8 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index_all
    3006              :       REAL(KIND=dp), DIMENSION(3), OPTIONAL              :: rcc
    3007              :       REAL(KIND=dp), DIMENSION(3)                        :: my_rcc
    3008              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
    3009            8 :       REAL(KIND=dp), DIMENSION(:), ALLOCATABLE           :: eigenvals
    3010            8 :       REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE        :: bc, xkp
    3011              :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
    3012              :          ALLOCATABLE, OPTIONAL                           :: berry_c
    3013              :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
    3014            8 :          ALLOCATABLE                                     :: D_rs, H_rs, S_rs
    3015              :       TYPE(cell_type), POINTER                           :: cell
    3016            8 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: moments_rs_img, matrix_ks_kp, &
    3017            8 :                                                             matrix_s_kp
    3018              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3019              :       TYPE(kpoint_type), POINTER                         :: kpoints_all, kpoints_scf
    3020            8 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    3021              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    3022              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    3023            8 :          POINTER                                         :: sab_all
    3024              : 
    3025            8 :       CALL timeset(routineN, handle)
    3026            8 :       calc_bc = PRESENT(berry_c)
    3027            8 :       my_do_parallel = .FALSE.
    3028            8 :       my_rcc = 0.0_dp
    3029            8 :       IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
    3030            8 :       IF (PRESENT(rcc)) my_rcc = rcc
    3031              : 
    3032              :       CALL get_qs_env(qs_env, &
    3033              :                       matrix_ks_kp=matrix_ks_kp, &
    3034              :                       matrix_s_kp=matrix_s_kp, &
    3035              :                       sab_all=sab_all, &
    3036              :                       cell=cell, &
    3037              :                       kpoints=kpoints_scf, &
    3038              :                       para_env=para_env, &
    3039              :                       dft_control=dft_control, &
    3040            8 :                       mos=mos)
    3041              : 
    3042            8 :       CALL get_mo_set(mo_set=mos(1), nao=nao)
    3043            8 :       CALL get_cell(cell=cell, h=hmat, periodic=periodic)
    3044            8 :       nspin = SIZE(matrix_ks_kp, 1)
    3045            8 :       nkp = SIZE(xkp, 2)
    3046              : 
    3047              :       ! create kpoint environment kpoints_all which contains all neighbor cells R
    3048              :       ! without considering any lattice symmetry
    3049            8 :       NULLIFY (kpoints_all)
    3050            8 :       CALL kpoint_create(kpoints_all)
    3051            8 :       CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_scf)
    3052              : 
    3053              :       CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, &
    3054            8 :                            index_to_cell=index_to_cell_all)
    3055            8 :       n_img_all = SIZE(index_to_cell_all, 2)
    3056              : 
    3057            8 :       NULLIFY (moments_rs_img)
    3058            8 :       CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_all)
    3059              :       ! D_μ,ν = <φ_μ|r|φ_ν>
    3060            8 :       CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
    3061              : 
    3062           80 :       ALLOCATE (S_rs(1, nao, nao, n_img_all), H_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
    3063           40 :       ALLOCATE (D_rs(3, nao, nao, n_img_all), source=0.0_dp)
    3064              : 
    3065              :       ! Convert real-space dbcsr matrices into arrays
    3066            8 :       CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, S_rs, cell_to_index_all)
    3067            8 :       CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, H_rs, cell_to_index_all)
    3068            8 :       CALL replicate_rs_matrices(moments_rs_img, kpoints_all, D_rs, cell_to_index_all)
    3069              : 
    3070            8 :       mepos = 0
    3071            8 :       num_pe = 1
    3072            8 :       num_copy = nkp
    3073            8 :       IF (my_do_parallel) THEN
    3074            8 :          mepos = para_env%mepos
    3075            8 :          num_pe = para_env%num_pe
    3076            8 :          num_copy = CEILING(REAL(nkp)/num_pe)
    3077              :       END IF
    3078              : 
    3079           56 :       ALLOCATE (dipole(nspin, num_copy, 3, nao, nao), source=z_zero)
    3080           32 :       IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp)
    3081              : 
    3082              : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, &
    3083              : !$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) &
    3084              : !$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, periodic, &
    3085            8 : !$OMP nkp, xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_bc)
    3086              :       ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero)
    3087              :       ALLOCATE (CDC(nao, nao), dip(3, nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero)
    3088              :       ALLOCATE (C_k(nao, nao), D_k(nao, nao), source=z_zero)
    3089              :       ALLOCATE (eigenvals(nao), source=0.0_dp)
    3090              :       IF (calc_bc) ALLOCATE (bc(3, nao), source=0.0_dp)
    3091              : !$OMP DO COLLAPSE(2)
    3092              :       DO ispin = 1, nspin
    3093              :          DO ikp = 1, nkp
    3094              :             IF (MOD(ikp - 1, num_pe) /= mepos) CYCLE
    3095              : 
    3096              :             ! S^R -> S(k), H^R -> H(k)
    3097              :             S_k = 0
    3098              :             H_k = 0
    3099              :             CALL rs_to_kp(S_rs(1, :, :, :), S_k, index_to_cell_all, xkp(:, ikp))
    3100              :             CALL rs_to_kp(H_rs(ispin, :, :, :), H_k, index_to_cell_all, xkp(:, ikp))
    3101              : 
    3102              :             ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
    3103              :             CALL geeig_right(H_k, S_k, eigenvals, C_k)
    3104              : 
    3105              :             ! To have a smooth complex phase of C(k) as function of k, for every n, we force
    3106              :             ! the largest C_μ,n(k) to be real.
    3107              :             ! This is important to have a continuous dipole moment d_nm(k) as a function of k
    3108              :             DO n = 1, nao
    3109              :                tmp_max = C_k(1, n)
    3110              :                DO mu = 1, nao
    3111              :                   IF (ABS(C_k(mu, n)) < ABS(tmp_max)) CYCLE
    3112              :                   tmp_max = C_k(mu, n)
    3113              :                END DO
    3114              :                phase = tmp_max/ABS(tmp_max)
    3115              :                C_k(:, n) = C_k(:, n)/phase
    3116              :             END DO
    3117              : 
    3118              :             DO i_dir = 1, 3 ! d^x, d^y, d^z
    3119              : 
    3120              :                IF (periodic(i_dir) == 0) CYCLE
    3121              :                ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR)
    3122              :                CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
    3123              :                CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
    3124              : 
    3125              :                ! Σ_R D^R e^(ikR) = D(k), D_μ,ν = <φ_μ|r|φ_ν>
    3126              :                CALL rs_to_kp(D_rs(i_dir, :, :, :), D_k(:, :), index_to_cell_all, xkp(:, ikp))
    3127              : 
    3128              :                ! Basis transform to Kohn-Sham basis: (C^H) ∇ S C, (C^H) ∇ H C, (C^H) D C
    3129              :                CALL gemm_square(C_k, 'C', dS_dk_i, 'N', C_k, 'N', C_dS_C)
    3130              :                CALL gemm_square(C_k, 'C', dH_dk_i, 'N', C_k, 'N', C_dH_C)
    3131              :                CALL gemm_square(C_k, 'C', D_k, 'N', C_k, 'N', CDC)
    3132              : 
    3133              :                ! Compute the dipole
    3134              :                ! d_nm (k) = - i/(ε(n)-ε(m)) [ (C^H)(dH(k)/dk)C ]_nm
    3135              :                !            + i ε(n)/(ε(n)-ε(m)) [ (C^H)(dS(k)/dk)C ]_nm + [ (C^H)D(k)C ]_nm
    3136              :                DO n = 1, nao
    3137              :                   DO m = 1, nao
    3138              :                      IF (n == m) CYCLE ! diagonal elements would need to be computed from
    3139              :                      ! a numerical k-derivative which is not implemented
    3140              :                      dip(i_dir, n, m) = -gaussi*C_dH_C(n, m)/(eigenvals(n) - eigenvals(m)) &
    3141              :                                         + gaussi*eigenvals(n)*C_dS_C(n, m)/(eigenvals(n) - eigenvals(m)) &
    3142              :                                         + CDC(n, m)
    3143              :                   END DO
    3144              :                END DO
    3145              :             END DO
    3146              :             ! Compute the Berry curvature from the dipoles
    3147              :             ! Ω^γ_n = Σ_m 2*Im[d^α_nm d^β_mn], where, α, β, γ belong to {x, y, z}
    3148              :             IF (calc_bc) THEN
    3149              :                bc = 0.0_dp
    3150              :                DO i_dir = 1, 3
    3151              :                   DO n = 1, nao
    3152              :                      DO m = 1, nao
    3153              :                         IF (n == m) CYCLE
    3154              :                         bc(i_dir, n) = bc(i_dir, n) &
    3155              :                                        + 2*AIMAG(dip(1 + MOD(i_dir, 3), n, m)*dip(1 + MOD(i_dir + 1, 3), m, n))
    3156              :                      END DO
    3157              :                   END DO
    3158              :                END DO
    3159              :             END IF
    3160              :             ! Store the dipoles and berry curvature for each MPI rank
    3161              :             dipole(ispin, CEILING(REAL(ikp)/num_pe), :, :, :) = dip(:, :, :)
    3162              :             IF (calc_bc) berry_c(ispin, CEILING(REAL(ikp)/num_pe), :, :) = bc(:, :)
    3163              :          END DO
    3164              :       END DO
    3165              : !$OMP END DO
    3166              :       DEALLOCATE (dS_dk_i, C_dS_C, dH_dk_i, C_dH_C, CDC, dip, S_k, H_k, C_k, D_k, eigenvals)
    3167              :       IF (calc_bc) DEALLOCATE (bc)
    3168              : !$OMP END PARALLEL
    3169            8 :       DEALLOCATE (S_rs, H_rs, D_rs)
    3170            8 :       CALL dbcsr_deallocate_matrix_set(moments_rs_img)
    3171            8 :       CALL kpoint_release(kpoints_all)
    3172            8 :       CALL timestop(handle)
    3173           32 :    END SUBROUTINE qs_moment_kpoints_deep
    3174              : 
    3175              : ! **************************************************************************************************
    3176              : !> \brief Calculates interband k-point dipoles in the existing SCF MO basis.
    3177              : !> \param qs_env ...
    3178              : !> \param dipole ...
    3179              : !> \param rcc retained for interface compatibility; interband dipoles are origin independent
    3180              : !> \param nmo_spin_out number of SCF MOs available for each spin
    3181              : ! **************************************************************************************************
    3182            6 :    SUBROUTINE qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_out)
    3183              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3184              :       COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
    3185              :          ALLOCATABLE                                     :: dipole
    3186              :       REAL(KIND=dp), DIMENSION(3), OPTIONAL              :: rcc
    3187              :       INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT), &
    3188              :          OPTIONAL                                        :: nmo_spin_out
    3189              : 
    3190              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_moment_kpoints_scf_mos'
    3191              : 
    3192              :       INTEGER                                            :: handle, i_dir, ikp, ikp_local, ispin, &
    3193              :                                                             m, n, nao, nkp, nmo, nspin
    3194            6 :       INTEGER, DIMENSION(:), ALLOCATABLE                 :: nmo_spin
    3195              :       INTEGER, DIMENSION(2)                              :: kp_range
    3196            6 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    3197              :       LOGICAL                                            :: my_kpgrp
    3198              :       REAL(KIND=dp), PARAMETER                           :: eps_degenerate = 1.0E-10_dp
    3199              :       REAL(KIND=dp)                                      :: cimag, creal, energy_diff
    3200            6 :       REAL(KIND=dp), DIMENSION(:), ALLOCATABLE           :: eigenvalues_kp
    3201            6 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvals
    3202              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_all
    3203              :       TYPE(cp_fm_struct_type), POINTER                   :: moment_struct
    3204              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    3205              :       TYPE(cp_fm_type)                                   :: fm_dummy, fm_tmp, mo_coeff_im_global, &
    3206              :                                                             mo_coeff_re_global, moment_im, &
    3207              :                                                             moment_re
    3208              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff_im, mo_coeff_re
    3209            6 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: overlap_deriv
    3210              :       TYPE(dbcsr_type), POINTER                          :: cmatrix, rmatrix
    3211              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3212            6 :       TYPE(kpoint_env_p_type), DIMENSION(:), POINTER     :: kp_env
    3213              :       TYPE(kpoint_env_type), POINTER                     :: kp
    3214              :       TYPE(kpoint_type), POINTER                         :: kpoints_scf
    3215            6 :       TYPE(mo_set_type), DIMENSION(:, :), POINTER        :: mos_kp
    3216              :       TYPE(mp_para_env_type), POINTER                    :: para_env, para_env_kp
    3217              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    3218            6 :          POINTER                                         :: sab_kp, sab_orb
    3219              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    3220              : 
    3221            6 :       CALL timeset(routineN, handle)
    3222              : 
    3223            6 :       NULLIFY (blacs_env_all, cell_to_index, cmatrix, dft_control, eigenvals, fm_struct, kp, &
    3224            6 :                kp_env, kpoints_scf, ks_env, mo_coeff_im, mo_coeff_re, moment_struct, mos_kp, &
    3225            6 :                overlap_deriv, para_env, para_env_kp, rmatrix, sab_kp, sab_orb)
    3226              :       IF (PRESENT(rcc)) THEN
    3227              :          MARK_USED(rcc)
    3228              :       END IF
    3229              : 
    3230              :       CALL get_qs_env(qs_env, dft_control=dft_control, kpoints=kpoints_scf, ks_env=ks_env, &
    3231            6 :                       para_env=para_env, sab_orb=sab_orb)
    3232            6 :       CPASSERT(ASSOCIATED(dft_control))
    3233            6 :       CPASSERT(ASSOCIATED(kpoints_scf))
    3234            6 :       CPASSERT(ASSOCIATED(ks_env))
    3235            6 :       CPASSERT(ASSOCIATED(para_env))
    3236            6 :       CPASSERT(ASSOCIATED(sab_orb))
    3237              : 
    3238              :       CALL get_kpoint_info(kpoints_scf, nkp=nkp, kp_range=kp_range, kp_env=kp_env, &
    3239              :                            para_env_kp=para_env_kp, blacs_env_all=blacs_env_all, &
    3240            6 :                            cell_to_index=cell_to_index, sab_nl=sab_kp)
    3241            6 :       IF (kp_range(2) >= kp_range(1)) THEN
    3242            6 :          CPASSERT(ASSOCIATED(kp_env))
    3243              :       END IF
    3244            6 :       CPASSERT(ASSOCIATED(para_env_kp))
    3245            6 :       CPASSERT(ASSOCIATED(blacs_env_all))
    3246            6 :       CPASSERT(ASSOCIATED(cell_to_index))
    3247            6 :       CPASSERT(ASSOCIATED(sab_kp))
    3248              : 
    3249              :       CALL build_overlap_matrix(ks_env, matrixkp_s=overlap_deriv, nderivative=1, &
    3250              :                                 basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
    3251            6 :                                 ext_kpoints=kpoints_scf)
    3252              : 
    3253            6 :       nspin = dft_control%nspins
    3254            6 :       CALL dbcsr_get_info(overlap_deriv(1, 1)%matrix, nfullrows_total=nao)
    3255           18 :       ALLOCATE (nmo_spin(nspin), source=0)
    3256            6 :       IF (kp_range(2) >= kp_range(1)) THEN
    3257            6 :          kp => kp_env(1)%kpoint_env
    3258            6 :          mos_kp => kp%mos
    3259            6 :          CPASSERT(ASSOCIATED(mos_kp))
    3260           12 :          DO ispin = 1, nspin
    3261           12 :             CALL get_mo_set(mos_kp(1, ispin), nmo=nmo_spin(ispin))
    3262              :          END DO
    3263              :       END IF
    3264            6 :       CALL para_env%max(nmo_spin)
    3265           48 :       ALLOCATE (dipole(nspin, nkp, 3, MAXVAL(nmo_spin), MAXVAL(nmo_spin)), source=z_zero)
    3266            6 :       IF (PRESENT(nmo_spin_out)) THEN
    3267            8 :          ALLOCATE (nmo_spin_out(nspin))
    3268            8 :          nmo_spin_out(:) = nmo_spin(:)
    3269              :       END IF
    3270              : 
    3271            6 :       ALLOCATE (rmatrix, cmatrix)
    3272              :       CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
    3273            6 :                         matrix_type=dbcsr_type_antisymmetric)
    3274              :       CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
    3275            6 :                         matrix_type=dbcsr_type_symmetric)
    3276            6 :       CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_kp)
    3277            6 :       CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_kp)
    3278              : 
    3279          258 :       DO ikp = 1, nkp
    3280          252 :          my_kpgrp = (ikp >= kp_range(1) .AND. ikp <= kp_range(2))
    3281              :          IF (my_kpgrp) THEN
    3282          234 :             ikp_local = ikp - kp_range(1) + 1
    3283          234 :             kp => kp_env(ikp_local)%kpoint_env
    3284          234 :             mos_kp => kp%mos
    3285              :          ELSE
    3286          252 :             NULLIFY (kp, mos_kp)
    3287              :          END IF
    3288          510 :          DO ispin = 1, nspin
    3289          252 :             nmo = nmo_spin(ispin)
    3290          756 :             ALLOCATE (eigenvalues_kp(nmo), source=0.0_dp)
    3291              : 
    3292              :             CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
    3293          252 :                                      para_env=para_env, context=blacs_env_all)
    3294          252 :             CALL cp_fm_create(mo_coeff_re_global, fm_struct)
    3295          252 :             CALL cp_fm_create(mo_coeff_im_global, fm_struct)
    3296          252 :             CALL cp_fm_create(fm_tmp, fm_struct)
    3297          252 :             CALL cp_fm_struct_release(fm_struct)
    3298              :             CALL cp_fm_struct_create(moment_struct, nrow_global=nmo, ncol_global=nmo, &
    3299          252 :                                      para_env=para_env, context=blacs_env_all)
    3300          252 :             CALL cp_fm_create(moment_re, moment_struct)
    3301          252 :             CALL cp_fm_create(moment_im, moment_struct)
    3302          252 :             CALL cp_fm_struct_release(moment_struct)
    3303              : 
    3304          252 :             IF (my_kpgrp) THEN
    3305          234 :                CALL get_mo_set(mos_kp(1, ispin), eigenvalues=eigenvals, mo_coeff=mo_coeff_re)
    3306          234 :                CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
    3307          234 :                CPASSERT(ASSOCIATED(eigenvals))
    3308          234 :                CPASSERT(ASSOCIATED(mo_coeff_re))
    3309          234 :                CPASSERT(ASSOCIATED(mo_coeff_im))
    3310          772 :                IF (para_env_kp%is_source()) eigenvalues_kp(1:nmo) = eigenvals(1:nmo)
    3311          234 :                CALL cp_fm_copy_general(mo_coeff_re, mo_coeff_re_global, para_env)
    3312          234 :                CALL cp_fm_copy_general(mo_coeff_im, mo_coeff_im_global, para_env)
    3313              :             ELSE
    3314           18 :                CALL cp_fm_copy_general(fm_dummy, mo_coeff_re_global, para_env)
    3315           18 :                CALL cp_fm_copy_general(fm_dummy, mo_coeff_im_global, para_env)
    3316              :             END IF
    3317          252 :             CALL para_env%sum(eigenvalues_kp)
    3318              : 
    3319         1008 :             DO i_dir = 1, 3
    3320          756 :                CALL dbcsr_set(rmatrix, 0.0_dp)
    3321          756 :                CALL dbcsr_set(cmatrix, 0.0_dp)
    3322              :                CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=overlap_deriv, &
    3323              :                                    ispin=i_dir + 1, xkp=kpoints_scf%xkp(:, ikp), &
    3324          756 :                                    cell_to_index=cell_to_index, sab_nl=sab_kp)
    3325              : 
    3326              :                ! Project the complex AO derivative operator as C^H A C. The
    3327              :                ! off-diagonal length-gauge dipoles follow from the energy-gap relation.
    3328          756 :                CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo)
    3329              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3330          756 :                                   1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
    3331              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3332          756 :                                   -1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)
    3333              : 
    3334          756 :                CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo)
    3335              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3336          756 :                                   1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
    3337              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3338          756 :                                   1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
    3339              : 
    3340          756 :                CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo)
    3341              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3342          756 :                                   1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
    3343              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3344          756 :                                   1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
    3345              : 
    3346          756 :                CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo)
    3347              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3348          756 :                                   -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
    3349              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
    3350          756 :                                   1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
    3351              : 
    3352         4236 :                DO n = 1, nmo
    3353        18108 :                   DO m = 1, nmo
    3354        14124 :                      IF (n == m) CYCLE
    3355        10896 :                      energy_diff = eigenvalues_kp(m) - eigenvalues_kp(n)
    3356        10896 :                      IF (ABS(energy_diff) <= eps_degenerate) CYCLE
    3357        10872 :                      CALL cp_fm_get_element(moment_re, m, n, creal)
    3358        10872 :                      CALL cp_fm_get_element(moment_im, m, n, cimag)
    3359        14100 :                      IF (para_env%is_source()) then
    3360         8688 :                         dipole(ispin, ikp, i_dir, n, m) = CMPLX(creal, cimag, KIND=dp)/energy_diff
    3361              :                      end if
    3362              :                   END DO
    3363              :                END DO
    3364              :             END DO
    3365          252 :             CALL cp_fm_release(mo_coeff_im_global)
    3366          252 :             CALL cp_fm_release(mo_coeff_re_global)
    3367          252 :             CALL cp_fm_release(moment_im)
    3368          252 :             CALL cp_fm_release(moment_re)
    3369          252 :             CALL cp_fm_release(fm_tmp)
    3370         1008 :             DEALLOCATE (eigenvalues_kp)
    3371              :          END DO
    3372              :       END DO
    3373              : 
    3374           12 :       DO ispin = 1, nspin
    3375          264 :          DO ikp = 1, nkp
    3376         1014 :             DO i_dir = 1, 3
    3377        35712 :                CALL para_env%sum(dipole(ispin, ikp, i_dir, :, :))
    3378              :             END DO
    3379              :          END DO
    3380              :       END DO
    3381              : 
    3382            6 :       CALL dbcsr_deallocate_matrix(cmatrix)
    3383            6 :       CALL dbcsr_deallocate_matrix(rmatrix)
    3384            6 :       CALL dbcsr_deallocate_matrix_set(overlap_deriv)
    3385            6 :       DEALLOCATE (nmo_spin)
    3386            6 :       CALL timestop(handle)
    3387              : 
    3388           18 :    END SUBROUTINE qs_moment_kpoints_scf_mos
    3389              : 
    3390              : ! **************************************************************************************************
    3391              : !> \brief Calculate and print dipole moment elements d_nm(k) for k-point calculations
    3392              : !> \param qs_env ...
    3393              : !> \param nmoments ...
    3394              : !> \param reference ...
    3395              : !> \param ref_point ...
    3396              : !> \param max_nmo ...
    3397              : !> \param unit_number ...
    3398              : !> \author Shridhar Shanbhag
    3399              : ! **************************************************************************************************
    3400           10 :    SUBROUTINE qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
    3401              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3402              :       INTEGER, INTENT(IN)                                :: nmoments, reference, max_nmo
    3403              :       REAL(dp), DIMENSION(:), INTENT(IN), POINTER        :: ref_point
    3404              :       INTEGER, INTENT(IN)                                :: unit_number
    3405              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_moment_kpoints'
    3406           10 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp
    3407           10 :       COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE  :: dipole_to_print
    3408              :       COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
    3409           10 :          ALLOCATABLE                                     :: dipole
    3410              :       INTEGER                                            :: handle, i_dir, ikp, nmo_dim, nkp, nao, &
    3411              :                                                             num_pe, mepos, n, m, &
    3412              :                                                             ispin, nspin, nmin, nmax, homo
    3413           10 :       INTEGER, DIMENSION(:), ALLOCATABLE                 :: nmo_spin_scf
    3414              :       LOGICAL                                            :: explicit_kpnts, explicit_kpset, use_scf_mos
    3415              :       REAL(KIND=dp), DIMENSION(3)                        :: rcc
    3416           10 :       REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE        :: xkp
    3417           10 :       REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE        :: bc_to_print
    3418           10 :       REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE  :: berry_c
    3419           10 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    3420              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    3421              :       TYPE(section_vals_type), POINTER                   :: kpnts, kpset
    3422              :       CHARACTER(LEN=default_string_length), &
    3423           10 :          DIMENSION(:), POINTER                           :: special_pnts
    3424              : 
    3425           10 :       CALL timeset(routineN, handle)
    3426              : 
    3427           10 :       IF (nmoments > 1) CPABORT("KPOINT quadrupole and higher moments not implemented.")
    3428           10 :       IF (max_nmo < 0) CPABORT("Negative maximum number of molecular orbitals max_nmo provided.")
    3429              : 
    3430              :       CALL get_qs_env(qs_env, &
    3431              :                       para_env=para_env, &
    3432              :                       matrix_ks_kp=matrix_ks_kp, &
    3433           10 :                       mos=mos)
    3434              : 
    3435           10 :       CALL get_mo_set(mo_set=mos(1), nao=nao)
    3436           10 :       CALL get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
    3437           10 :       CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
    3438           10 :       nspin = SIZE(matrix_ks_kp, 1)
    3439           10 :       nkp = SIZE(xkp, 2)
    3440              : 
    3441           10 :       kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
    3442           10 :       kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
    3443           10 :       CALL section_vals_get(kpset, explicit=explicit_kpset)
    3444           10 :       CALL section_vals_get(kpnts, explicit=explicit_kpnts)
    3445           10 :       use_scf_mos = .NOT. explicit_kpset .AND. .NOT. explicit_kpnts
    3446              : 
    3447           10 :       IF (unit_number > 0) WRITE (unit_number, FMT="(/,T2,A)") &
    3448            5 :          '!-----------------------------------------------------------------------------!'
    3449           10 :       IF (unit_number > 0) WRITE (unit_number, "(T22,A)") "Periodic Dipole Matrix Elements"
    3450              : 
    3451           10 :       IF (use_scf_mos) THEN
    3452            4 :          CALL qs_moment_kpoints_scf_mos(qs_env, dipole, rcc, nmo_spin_scf)
    3453            4 :          nmo_dim = SIZE(dipole, 4)
    3454           24 :          ALLOCATE (berry_c(nspin, nkp, 3, nmo_dim), source=0.0_dp)
    3455            8 :          DO ispin = 1, nspin
    3456          256 :             DO ikp = 1, nkp
    3457          996 :                DO i_dir = 1, 3
    3458         4160 :                   DO n = 1, nmo_dim
    3459        17736 :                      DO m = 1, nmo_dim
    3460        13824 :                         IF (n == m) CYCLE
    3461              :                         berry_c(ispin, ikp, i_dir, n) = berry_c(ispin, ikp, i_dir, n) &
    3462              :                                                         + 2*AIMAG(dipole(ispin, ikp, 1 + MOD(i_dir, 3), n, m)* &
    3463        16992 :                                                                   dipole(ispin, ikp, 1 + MOD(i_dir + 1, 3), m, n))
    3464              :                      END DO
    3465              :                   END DO
    3466              :                END DO
    3467              :             END DO
    3468              :          END DO
    3469              :       ELSE
    3470              :          CALL qs_moment_kpoints_deep(qs_env, &
    3471              :                                      xkp, &
    3472              :                                      dipole, &
    3473              :                                      rcc, &
    3474              :                                      berry_c, &
    3475            6 :                                      do_parallel=.TRUE.)
    3476            6 :          nmo_dim = nao
    3477              :       END IF
    3478              : 
    3479           40 :       ALLOCATE (dipole_to_print(3, nmo_dim, nmo_dim), source=z_zero)
    3480           30 :       ALLOCATE (bc_to_print(3, nmo_dim), source=0.0_dp)
    3481              : 
    3482           10 :       mepos = para_env%mepos
    3483           10 :       num_pe = para_env%num_pe
    3484              : 
    3485          264 :       DO ikp = 1, nkp
    3486          520 :          DO ispin = 1, nspin
    3487          256 :             CALL get_mo_set(mo_set=mos(ispin), homo=homo)
    3488          256 :             nmin = max(1, homo - (max_nmo - 1)/2)
    3489          256 :             nmax = min(nao, homo + max_nmo/2)
    3490          256 :             IF (max_nmo == 0) THEN
    3491            0 :                nmin = 1
    3492            0 :                nmax = nao
    3493              :             END IF
    3494          256 :             IF (use_scf_mos) THEN
    3495          248 :                nmax = min(nmax, nmo_spin_scf(ispin))
    3496              :             END IF
    3497          256 :             dipole_to_print = 0.0_dp
    3498          256 :             bc_to_print = 0.0_dp
    3499          256 :             IF (use_scf_mos) THEN
    3500        19736 :                dipole_to_print(:, :, :) = dipole(ispin, ikp, :, :, :)
    3501         4472 :                bc_to_print(:, :) = berry_c(ispin, ikp, :, :)
    3502            8 :             ELSE IF (mod(ikp - 1, num_pe) == mepos) THEN
    3503        87268 :                dipole_to_print(:, :, :) = dipole(ispin, CEILING(REAL(ikp)/num_pe), :, :, :)
    3504          900 :                bc_to_print(:, :) = berry_c(ispin, CEILING(REAL(ikp)/num_pe), :, :)
    3505              :             END IF
    3506          256 :             IF (.NOT. use_scf_mos) THEN
    3507            8 :                CALL para_env%sum(dipole_to_print)
    3508            8 :                CALL para_env%sum(bc_to_print)
    3509              :             END IF
    3510          766 :             IF (unit_number > 0) THEN
    3511          128 :                IF (special_pnts(ikp) /= "") WRITE (unit_number, "(/,2X,A,A)") &
    3512            3 :                   "Special point: ", ADJUSTL(TRIM(special_pnts(ikp)))
    3513              :                WRITE (unit_number, "(/,1X,A,I3,1X,3(A,1F12.6))") &
    3514          128 :                   "Kpoint:", ikp, ", kx:", xkp(1, ikp), ", ky:", xkp(2, ikp), ", kz:", xkp(3, ikp)
    3515          128 :                IF (nspin > 1) WRITE (unit_number, "(/,2X,A,I2)") "Open Shell System. Spin:", ispin
    3516              :                WRITE (unit_number, "(2X,A)") "  kp   n   m  Re(dx_nm)  Im(dx_nm) &
    3517          128 :                & Re(dy_nm)  Im(dy_nm)  Re(dz_nm)  Im(dz_nm)"
    3518          676 :                DO n = nmin, nmax
    3519         3116 :                   DO m = nmin, nmax
    3520         2440 :                      IF (n == m) CYCLE
    3521         2988 :                      WRITE (unit_number, "(2X,I4,2I4,6(G11.3))") ikp, n, m, dipole_to_print(1:3, n, m)
    3522              :                   END DO
    3523              :                END DO
    3524          128 :                WRITE (unit_number, "(/,1X,A)") "Berry Curvature"
    3525          128 :                WRITE (unit_number, "(2X,A)") "   kp    n      YZ          ZX          XY"
    3526          676 :                DO n = nmin, nmax
    3527              :                   WRITE (unit_number, "(2X,2I5,3(1X,G11.3))") &
    3528          676 :                      ikp, n, bc_to_print(1, n), bc_to_print(2, n), bc_to_print(3, n)
    3529              :                END DO
    3530              :             END IF
    3531              :          END DO
    3532              :       END DO
    3533           10 :       DEALLOCATE (dipole_to_print, bc_to_print, berry_c, dipole)
    3534           10 :       IF (ALLOCATED(nmo_spin_scf)) DEALLOCATE (nmo_spin_scf)
    3535           10 :       DEALLOCATE (special_pnts, xkp)
    3536              : 
    3537           10 :       CALL timestop(handle)
    3538              : 
    3539           40 :    END SUBROUTINE qs_moment_kpoints
    3540              : 
    3541              : END MODULE qs_moments
        

Generated by: LCOV version 2.0-1