LCOV - code coverage report
Current view: top level - src - qs_operators_ao.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.0 % 682 675
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

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

Generated by: LCOV version 2.0-1