LCOV - code coverage report
Current view: top level - src - qs_vcd_ao.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 93.8 % 536 503
Test Date: 2025-07-25 12:55:17 Functions: 92.3 % 13 12

            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              : MODULE qs_vcd_ao
       8              :    USE ai_contraction,                  ONLY: block_add,&
       9              :                                               contraction
      10              :    USE ai_kinetic,                      ONLY: kinetic
      11              :    USE ai_overlap_ppl,                  ONLY: ppl_integral
      12              :    USE ao_util,                         ONLY: exp_radius_very_extended
      13              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14              :                                               get_atomic_kind_set
      15              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      16              :                                               gto_basis_set_p_type,&
      17              :                                               gto_basis_set_type
      18              :    USE block_p_types,                   ONLY: block_p_type
      19              :    USE cell_types,                      ONLY: cell_type,&
      20              :                                               pbc
      21              :    USE cp_control_types,                ONLY: dft_control_type
      22              :    USE cp_dbcsr_api,                    ONLY: &
      23              :         dbcsr_add, dbcsr_copy, dbcsr_desymmetrize, dbcsr_distribution_get, &
      24              :         dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_init_p, &
      25              :         dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_work_create
      26              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      27              :                                               dbcsr_deallocate_matrix_set
      28              :    USE external_potential_types,        ONLY: get_potential,&
      29              :                                               gth_potential_type,&
      30              :                                               sgp_potential_type
      31              :    USE gaussian_gridlevels,             ONLY: gridlevel_info_type
      32              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      33              :                                               section_vals_type
      34              :    USE kinds,                           ONLY: default_string_length,&
      35              :                                               dp,&
      36              :                                               int_8
      37              :    USE memory_utilities,                ONLY: reallocate
      38              :    USE message_passing,                 ONLY: mp_comm_type
      39              :    USE orbital_pointers,                ONLY: coset,&
      40              :                                               init_orbital_pointers,&
      41              :                                               ncoset
      42              :    USE particle_types,                  ONLY: particle_type
      43              :    USE pw_env_types,                    ONLY: pw_env_get,&
      44              :                                               pw_env_type
      45              :    USE pw_methods,                      ONLY: pw_axpy,&
      46              :                                               pw_zero
      47              :    USE pw_pool_types,                   ONLY: pw_pool_type
      48              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      49              :    USE qs_energy_types,                 ONLY: qs_energy_type
      50              :    USE qs_environment_types,            ONLY: get_qs_env,&
      51              :                                               qs_environment_type
      52              :    USE qs_integral_utils,               ONLY: basis_set_list_setup,&
      53              :                                               get_memory_usage
      54              :    USE qs_integrate_potential,          ONLY: integrate_pgf_product
      55              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      56              :                                               get_qs_kind_set,&
      57              :                                               qs_kind_type
      58              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      59              :    USE qs_linres_types,                 ONLY: vcd_env_type
      60              :    USE qs_moments,                      ONLY: build_dsdv_moments
      61              :    USE qs_neighbor_list_types,          ONLY: &
      62              :         get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
      63              :         neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
      64              :         neighbor_list_iterator_release, neighbor_list_set_p_type, nl_set_sub_iterator, &
      65              :         nl_sub_iterate
      66              :    USE qs_rho_types,                    ONLY: qs_rho_type
      67              :    USE qs_vxc,                          ONLY: qs_vxc_create
      68              :    USE realspace_grid_types,            ONLY: realspace_grid_type
      69              :    USE rs_pw_interface,                 ONLY: potential_pw2rs
      70              :    USE sap_kind_types,                  ONLY: alist_type,&
      71              :                                               build_sap_ints,&
      72              :                                               get_alist,&
      73              :                                               release_sap_int,&
      74              :                                               sap_int_type,&
      75              :                                               sap_sort
      76              :    USE task_list_types,                 ONLY: atom_pair_type,&
      77              :                                               task_list_type,&
      78              :                                               task_type
      79              :    USE virial_types,                    ONLY: virial_type
      80              : 
      81              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      82              : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
      83              : !$                    omp_init_lock, omp_set_lock, &
      84              : !$                    omp_unset_lock, omp_destroy_lock
      85              : 
      86              : #include "./base/base_uses.f90"
      87              : 
      88              :    IMPLICIT NONE
      89              : 
      90              :    PRIVATE
      91              : 
      92              : ! *** Global parameters ***
      93              : 
      94              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd_ao'
      95              :    INTEGER, PARAMETER :: bi_1 = 1, bi_x = 2, bi_y = 3, bi_z = 4, bi_xx = 5, &
      96              :                          bi_xy = 6, bi_xz = 7, bi_yy = 8, bi_yz = 9, bi_zz = 10
      97              : 
      98              : ! *** Public subroutines ***
      99              : 
     100              :    PUBLIC :: build_dSdV_matrix, build_com_rpnl_r, &
     101              :              hr_mult_by_delta_3d, hr_mult_by_delta_1d, build_dcom_rpnl, build_drpnl_matrix, &
     102              :              build_tr_matrix, &
     103              :              build_rcore_matrix, &
     104              :              build_rpnl_matrix, &
     105              :              build_matrix_r_vhxc, &
     106              :              build_matrix_hr_rh
     107              : 
     108              : CONTAINS
     109              : 
     110              : ! **************************************************************************************************
     111              : !> \brief Build the matrix Hr*delta_nu^\lambda - rH*delta_mu^\lambda
     112              : !> \param vcd_env ...
     113              : !> \param qs_env ...
     114              : !> \param rc ...
     115              : !> \author Edward Ditler
     116              : ! **************************************************************************************************
     117           36 :    SUBROUTINE build_matrix_hr_rh(vcd_env, qs_env, rc)
     118              :       TYPE(vcd_env_type)                                 :: vcd_env
     119              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     120              :       REAL(dp), DIMENSION(3)                             :: rc
     121              : 
     122              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_matrix_hr_rh'
     123              :       INTEGER, PARAMETER                                 :: ispin = 1
     124              : 
     125              :       INTEGER                                            :: handle, i
     126              :       REAL(KIND=dp)                                      :: eps_ppnl
     127              :       TYPE(cell_type), POINTER                           :: cell
     128              :       TYPE(dft_control_type), POINTER                    :: dft_control
     129              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     130           36 :          POINTER                                         :: sab_all, sap_ppnl
     131           36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     132           36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     133              : 
     134           36 :       CALL timeset(routineN, handle)
     135              : 
     136              :       CALL get_qs_env(qs_env=qs_env, &
     137              :                       dft_control=dft_control, &
     138              :                       particle_set=particle_set, &
     139              :                       sab_all=sab_all, &
     140              :                       sap_ppnl=sap_ppnl, &
     141              :                       qs_kind_set=qs_kind_set, &
     142           36 :                       cell=cell)
     143              : 
     144           36 :       eps_ppnl = dft_control%qs_control%eps_ppnl
     145          144 :       DO i = 1, 3
     146          108 :          CALL dbcsr_set(vcd_env%matrix_hr(ispin, i)%matrix, 0._dp)
     147          144 :          CALL dbcsr_set(vcd_env%matrix_rh(ispin, i)%matrix, 0._dp)
     148              :       END DO
     149              : 
     150              :       ASSOCIATE (my_matrix_hr_1d => vcd_env%matrix_hr(ispin, 1:3))
     151              :          CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
     152              :                                 dft_control%qs_control%eps_ppnl, cell, rc, &
     153           36 :                                 direction_Or=.TRUE.)
     154              :          CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
     155           36 :                               direction_Or=.TRUE., rc=rc)
     156           36 :          CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, rc)
     157           72 :          CALL build_matrix_r_vhxc(vcd_env%matrix_hr, qs_env, rc)
     158              :       END ASSOCIATE
     159              : 
     160              :       ASSOCIATE (my_matrix_hr_1d => vcd_env%matrix_rh(ispin, 1:3))
     161              :          CALL build_rpnl_matrix(my_matrix_hr_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
     162              :                                 dft_control%qs_control%eps_ppnl, cell, rc, &
     163           36 :                                 direction_Or=.FALSE.)
     164              :          CALL build_tr_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, &
     165           36 :                               direction_Or=.FALSE., rc=rc)
     166           36 :          CALL build_rcore_matrix(my_matrix_hr_1d, qs_env, qs_kind_set, "ORB", sab_all, rc)
     167           72 :          CALL build_matrix_r_vhxc(vcd_env%matrix_rh, qs_env, rc)
     168              :       END ASSOCIATE
     169              : 
     170           36 :       CALL timestop(handle)
     171           36 :    END SUBROUTINE build_matrix_hr_rh
     172              : 
     173              : ! **************************************************************************************************
     174              : !> \brief Product of r with V_nl. Adapted from build_com_rpnl.
     175              : !> \param matrix_rv ...
     176              : !> \param qs_kind_set ...
     177              : !> \param particle_set ...
     178              : !> \param sab_all ...
     179              : !> \param sap_ppnl ...
     180              : !> \param eps_ppnl ...
     181              : !> \param cell ...
     182              : !> \param ref_point ...
     183              : !> \param direction_Or ...
     184              : !> \author Edward Ditler, Tomas Zimmermann
     185              : ! **************************************************************************************************
     186           76 :    SUBROUTINE build_rpnl_matrix(matrix_rv, qs_kind_set, particle_set, sab_all, sap_ppnl, eps_ppnl, &
     187              :                                 cell, ref_point, direction_Or)
     188              : 
     189              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: matrix_rv
     190              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     191              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     192              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     193              :          POINTER                                         :: sab_all, sap_ppnl
     194              :       REAL(KIND=dp), INTENT(IN)                          :: eps_ppnl
     195              :       TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER     :: cell
     196              :       REAL(KIND=dp), DIMENSION(3)                        :: ref_point
     197              :       LOGICAL                                            :: direction_Or
     198              : 
     199              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_rpnl_matrix'
     200              : 
     201              :       INTEGER                                            :: handle, i, iab, iac, iatom, ibc, icol, &
     202              :                                                             ikind, irow, jatom, jkind, kac, kbc, &
     203              :                                                             kkind, na, natom, nb, nkind, np, slot
     204              :       INTEGER, DIMENSION(3)                              :: cell_b
     205              :       LOGICAL                                            :: found, ppnl_present
     206           76 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: achint, acint, bchint, bcint
     207              :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
     208          304 :       TYPE(block_p_type), DIMENSION(3)                   :: blocks_rv
     209           76 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
     210              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     211           76 :       TYPE(sap_int_type), DIMENSION(:), POINTER          :: sap_int
     212              : 
     213              : !$    INTEGER(kind=omp_lock_kind), &
     214           76 : !$       ALLOCATABLE, DIMENSION(:) :: locks
     215              : !$    INTEGER                                            :: lock_num, hash
     216              : !$    INTEGER, PARAMETER                                 :: nlock = 501
     217              : 
     218           76 :       ppnl_present = ASSOCIATED(sap_ppnl)
     219           76 :       IF (.NOT. ppnl_present) RETURN
     220              : 
     221           76 :       CALL timeset(routineN, handle)
     222           76 :       nkind = SIZE(qs_kind_set)
     223           76 :       natom = SIZE(particle_set)
     224              : 
     225              :       ! sap_int needs to be shared as multiple threads need to access this
     226           76 :       NULLIFY (sap_int)
     227          532 :       ALLOCATE (sap_int(nkind*nkind))
     228          380 :       DO i = 1, nkind*nkind
     229          304 :          NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
     230          380 :          sap_int(i)%nalist = 0
     231              :       END DO
     232              : 
     233              :       MARK_USED(ref_point)
     234              :       ! "nder" in moment_mode is "order"
     235              :       CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=1, moment_mode=.TRUE., &
     236           76 :                           particle_set=particle_set, cell=cell, refpoint=ref_point)
     237              : 
     238              :       ! *** Set up a sorting index
     239           76 :       CALL sap_sort(sap_int)
     240              : 
     241          380 :       ALLOCATE (basis_set(nkind))
     242          228 :       DO ikind = 1, nkind
     243          152 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     244          228 :          IF (ASSOCIATED(orb_basis_set)) THEN
     245          152 :             basis_set(ikind)%gto_basis_set => orb_basis_set
     246              :          ELSE
     247            0 :             NULLIFY (basis_set(ikind)%gto_basis_set)
     248              :          END IF
     249              :       END DO
     250              : 
     251              :       ! *** All integrals needed have been calculated and stored in sap_int
     252              :       ! *** We now calculate the commutator matrix elements
     253              : 
     254              : !$OMP PARALLEL &
     255              : !$OMP DEFAULT (NONE) &
     256              : !$OMP SHARED  (basis_set, matrix_rv, &
     257              : !$OMP          sap_int, nkind, eps_ppnl, locks, sab_all, direction_Or) &
     258              : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, &
     259              : !$OMP          iab, irow, icol, blocks_rv, &
     260              : !$OMP          found, iac, ibc, alist_ac, alist_bc, &
     261              : !$OMP          na, np, nb, kkind, kac, kbc, i, lock_num, &
     262           76 : !$OMP          hash, natom, acint, bcint, achint, bchint)
     263              : 
     264              : !$OMP SINGLE
     265              : !$    ALLOCATE (locks(nlock))
     266              : !$OMP END SINGLE
     267              : 
     268              : !$OMP DO
     269              : !$    DO lock_num = 1, nlock
     270              : !$       call omp_init_lock(locks(lock_num))
     271              : !$    END DO
     272              : !$OMP END DO
     273              : 
     274              : !$OMP DO SCHEDULE(GUIDED)
     275              : 
     276              :       DO slot = 1, sab_all(1)%nl_size
     277              : 
     278              :          ikind = sab_all(1)%nlist_task(slot)%ikind
     279              :          jkind = sab_all(1)%nlist_task(slot)%jkind
     280              :          iatom = sab_all(1)%nlist_task(slot)%iatom
     281              :          jatom = sab_all(1)%nlist_task(slot)%jatom
     282              :          cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
     283              : 
     284              :          IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
     285              :          IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
     286              :          iab = ikind + nkind*(jkind - 1)
     287              : 
     288              :          ! *** Create matrix blocks for a new matrix block column ***
     289              :          irow = iatom
     290              :          icol = jatom
     291              :          DO i = 1, 3
     292              :             CALL dbcsr_get_block_p(matrix_rv(i)%matrix, irow, icol, blocks_rv(i)%block, found)
     293              :             CPASSERT(found)
     294              :          END DO
     295              : 
     296              :          ! loop over all kinds for projector atom
     297              :          DO kkind = 1, nkind
     298              :             iac = ikind + nkind*(kkind - 1)
     299              :             ibc = jkind + nkind*(kkind - 1)
     300              :             IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
     301              :             IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
     302              :             CALL get_alist(sap_int(iac), alist_ac, iatom)
     303              :             CALL get_alist(sap_int(ibc), alist_bc, jatom)
     304              : 
     305              :             IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
     306              :             IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
     307              :             DO kac = 1, alist_ac%nclist
     308              :                DO kbc = 1, alist_bc%nclist
     309              :                   IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
     310              :                   IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
     311              :                      IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
     312              :                      acint => alist_ac%clist(kac)%acint
     313              :                      bcint => alist_bc%clist(kbc)%acint
     314              :                      achint => alist_ac%clist(kac)%achint
     315              :                      bchint => alist_bc%clist(kbc)%achint
     316              :                      na = SIZE(acint, 1)
     317              :                      np = SIZE(acint, 2)
     318              :                      nb = SIZE(bcint, 1)
     319              : !$                   hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
     320              : !$                   CALL omp_set_lock(locks(hash))
     321              :                      IF (direction_Or) THEN
     322              :                         ! Vnl*r
     323              :                         blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
     324              :                                                          MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2)))
     325              :                         blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
     326              :                                                          MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3)))
     327              :                         blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
     328              :                                                          MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4)))
     329              :                      ELSE
     330              :                         ! r*Vnl
     331              :                         blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
     332              :                                                          MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))
     333              :                         blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
     334              :                                                          MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))
     335              :                         blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
     336              :                                                          MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1)))
     337              :                      END IF
     338              : !$                   CALL omp_unset_lock(locks(hash))
     339              :                      EXIT ! We have found a match and there can be only one single match
     340              :                   END IF
     341              :                END DO
     342              :             END DO
     343              :          END DO
     344              :          DO i = 1, 3
     345              :             NULLIFY (blocks_rv(i)%block)
     346              :          END DO
     347              :       END DO
     348              : 
     349              : !$OMP DO
     350              : !$    DO lock_num = 1, nlock
     351              : !$       call omp_destroy_lock(locks(lock_num))
     352              : !$    END DO
     353              : !$OMP END DO
     354              : 
     355              : !$OMP SINGLE
     356              : !$    DEALLOCATE (locks)
     357              : !$OMP END SINGLE NOWAIT
     358              : 
     359              : !$OMP END PARALLEL
     360              : 
     361           76 :       CALL release_sap_int(sap_int)
     362              : 
     363           76 :       DEALLOCATE (basis_set)
     364              : 
     365           76 :       CALL timestop(handle)
     366              : 
     367          228 :    END SUBROUTINE build_rpnl_matrix
     368              : 
     369              : ! **************************************************************************************************
     370              : !> \brief   Calculation of the product Tr or rT over Cartesian Gaussian functions.
     371              : !> \param   matrix_tr ...
     372              : !> \param   qs_env ...
     373              : !> \param   qs_kind_set ...
     374              : !> \param   basis_type basis set to be used
     375              : !> \param   sab_nl pair list (must be consistent with basis sets!)
     376              : !> \param   direction_Or ...
     377              : !> \param   rc ...
     378              : !> \date    11.10.2010
     379              : !> \par     History
     380              : !>          Ported from qs_overlap, replaces code in build_core_hamiltonian
     381              : !>          Refactoring [07.2014] JGH
     382              : !>          Simplify options and use new kinetic energy integral routine
     383              : !>          Adapted from qs_kinetic [07.2016]
     384              : !>          Adapted from build_com_tr_matrix [2021] by ED
     385              : !> \author  JGH
     386              : !> \version 1.0
     387              : ! **************************************************************************************************
     388           76 :    SUBROUTINE build_tr_matrix(matrix_tr, qs_env, qs_kind_set, basis_type, sab_nl, direction_Or, rc)
     389              : 
     390              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: matrix_tr
     391              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     392              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     393              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     394              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     395              :          POINTER                                         :: sab_nl
     396              :       LOGICAL                                            :: direction_Or
     397              :       REAL(KIND=dp), DIMENSION(3), OPTIONAL              :: rc
     398              : 
     399              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'build_tr_matrix'
     400              : 
     401              :       INTEGER                                            :: handle, iatom, icol, ikind, ir, irow, &
     402              :                                                             iset, jatom, jkind, jset, ldsab, ltab, &
     403              :                                                             natom, ncoa, ncob, nkind, nseta, &
     404              :                                                             nsetb, sgfa, sgfb, slot
     405              :       INTEGER, DIMENSION(3)                              :: cell
     406           76 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     407           76 :                                                             npgfb, nsgfa, nsgfb
     408           76 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     409              :       LOGICAL                                            :: do_symmetric, found, trans
     410              :       REAL(KIND=dp)                                      :: tab
     411           76 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: qab, tkab
     412           76 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: kab
     413              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, rb, rbc
     414           76 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     415           76 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: kx_block, ky_block, kz_block, rpgfa, &
     416           76 :                                                             rpgfb, scon_a, scon_b, zeta, zetb
     417              :       TYPE(cell_type), POINTER                           :: qs_cell
     418           76 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     419              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     420           76 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     421              : 
     422              : !$    INTEGER(kind=omp_lock_kind), &
     423           76 : !$       ALLOCATABLE, DIMENSION(:) :: locks
     424              : !$    INTEGER                                            :: lock_num, hash, hash1, hash2
     425              : !$    INTEGER(KIND=int_8)                                :: iatom8
     426              : !$    INTEGER, PARAMETER                                 :: nlock = 501
     427              : 
     428              :       MARK_USED(int_8)
     429              : 
     430           76 :       CALL timeset(routineN, handle)
     431              : 
     432           76 :       nkind = SIZE(qs_kind_set)
     433              : 
     434              :       ! check for symmetry
     435           76 :       CPASSERT(SIZE(sab_nl) > 0)
     436           76 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
     437              : 
     438              :       ! prepare basis set
     439          380 :       ALLOCATE (basis_set_list(nkind))
     440           76 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
     441              : 
     442              :       ! *** Allocate work storage ***
     443           76 :       ldsab = get_memory_usage(qs_kind_set, basis_type)
     444              : 
     445              :       CALL get_qs_env(qs_env=qs_env, &
     446              :                       particle_set=particle_set, &
     447              :                       cell=qs_cell, &
     448           76 :                       natom=natom)
     449              : 
     450              : !$OMP PARALLEL DEFAULT(NONE) &
     451              : !$OMP SHARED (ldsab,do_symmetric, sab_nl, rc,&
     452              : !$OMP         ncoset,matrix_tr,basis_set_list,qs_cell,natom,locks)  &
     453              : !$OMP PRIVATE (kx_block,ky_block,kz_block,kab,qab,tab,ikind,jkind,iatom,jatom,rab,rac,rbc,cell, &
     454              : !$OMP          basis_set_a, basis_set_b, nseta, ncoa, ncob, ltab, nsetb, tkab, &
     455              : !$OMP          irow, icol, found, trans, sgfa, sgfb, iset, jset, &
     456              : !$OMP          hash, hash1, hash2, iatom8, slot, lock_num) &
     457              : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, zeta, scon_a) &
     458              : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, zetb, scon_b) &
     459              : !$OMP SHARED(particle_set, direction_or) &
     460           76 : !$OMP PRIVATE(ra, rb)
     461              : 
     462              : !$OMP SINGLE
     463              : !$    ALLOCATE (locks(nlock))
     464              : !$OMP END SINGLE
     465              : 
     466              : !$OMP DO
     467              : !$    DO lock_num = 1, nlock
     468              : !$       call omp_init_lock(locks(lock_num))
     469              : !$    END DO
     470              : !$OMP END DO
     471              : 
     472              :       ALLOCATE (kab(ldsab, ldsab, 3), qab(ldsab, ldsab))
     473              : 
     474              : !$OMP DO SCHEDULE(GUIDED)
     475              :       DO slot = 1, sab_nl(1)%nl_size
     476              : 
     477              :          ikind = sab_nl(1)%nlist_task(slot)%ikind
     478              :          jkind = sab_nl(1)%nlist_task(slot)%jkind
     479              :          iatom = sab_nl(1)%nlist_task(slot)%iatom
     480              :          jatom = sab_nl(1)%nlist_task(slot)%jatom
     481              :          cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
     482              :          rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
     483              : 
     484              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     485              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     486              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     487              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     488              : 
     489              : !$       iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
     490              : !$       hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
     491              : 
     492              :          ! basis ikind
     493              :          first_sgfa => basis_set_a%first_sgf
     494              :          la_max => basis_set_a%lmax
     495              :          la_min => basis_set_a%lmin
     496              :          npgfa => basis_set_a%npgf
     497              :          nsgfa => basis_set_a%nsgf_set
     498              :          rpgfa => basis_set_a%pgf_radius
     499              :          set_radius_a => basis_set_a%set_radius
     500              :          scon_a => basis_set_a%scon
     501              :          zeta => basis_set_a%zet
     502              :          ! basis jkind
     503              :          first_sgfb => basis_set_b%first_sgf
     504              :          lb_max => basis_set_b%lmax
     505              :          lb_min => basis_set_b%lmin
     506              :          npgfb => basis_set_b%npgf
     507              :          nsgfb => basis_set_b%nsgf_set
     508              :          rpgfb => basis_set_b%pgf_radius
     509              :          set_radius_b => basis_set_b%set_radius
     510              :          scon_b => basis_set_b%scon
     511              :          zetb => basis_set_b%zet
     512              : 
     513              :          nseta = basis_set_a%nset
     514              :          nsetb = basis_set_b%nset
     515              : 
     516              :          IF (do_symmetric) THEN
     517              :             IF (iatom <= jatom) THEN
     518              :                irow = iatom
     519              :                icol = jatom
     520              :             ELSE
     521              :                irow = jatom
     522              :                icol = iatom
     523              :             END IF
     524              :          ELSE
     525              :             irow = iatom
     526              :             icol = jatom
     527              :          END IF
     528              :          NULLIFY (kx_block)
     529              :          CALL dbcsr_get_block_p(matrix=matrix_tr(1)%matrix, &
     530              :                                 row=irow, col=icol, BLOCK=kx_block, found=found)
     531              :          CPASSERT(found)
     532              :          NULLIFY (ky_block)
     533              :          CALL dbcsr_get_block_p(matrix=matrix_tr(2)%matrix, &
     534              :                                 row=irow, col=icol, BLOCK=ky_block, found=found)
     535              :          CPASSERT(found)
     536              :          NULLIFY (kz_block)
     537              :          CALL dbcsr_get_block_p(matrix=matrix_tr(3)%matrix, &
     538              :                                 row=irow, col=icol, BLOCK=kz_block, found=found)
     539              :          CPASSERT(found)
     540              : 
     541              :          ! The kinetic integrals depend only on rab (also for the screening)
     542              :          tab = NORM2(rab)
     543              : 
     544              :          ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
     545              :          ra = pbc(particle_set(iatom)%r(:), qs_cell)
     546              :          rb(:) = ra(:) + rab(:)
     547              :          rac = pbc(rc, ra, qs_cell)
     548              :          rbc = rac + rab
     549              : 
     550              :          trans = do_symmetric .AND. (iatom > jatom)
     551              : 
     552              :          DO iset = 1, nseta
     553              : 
     554              :             ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     555              :             sgfa = first_sgfa(1, iset)
     556              : 
     557              :             DO jset = 1, nsetb
     558              : 
     559              :                IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
     560              : 
     561              : !$             hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
     562              : !$             hash = MOD(hash1 + hash2, nlock) + 1
     563              : 
     564              :                ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     565              :                sgfb = first_sgfb(1, jset)
     566              : 
     567              :                ! calculate integrals
     568              :                ltab = MAX(npgfa(iset)*ncoset(la_max(iset) + 1), npgfb(jset)*ncoset(lb_max(jset) + 1))
     569              :                ALLOCATE (tkab(ltab, ltab))
     570              :                CALL kinetic(la_max(iset) + 1, la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     571              :                             lb_max(jset) + 1, lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     572              :                             rab, tkab)
     573              :                ! commutator
     574              :                CALL ab_opr(la_max(iset), npgfa(iset), rpgfa(:, iset), la_min(iset), &
     575              :                            lb_max(jset), npgfb(jset), rpgfb(:, jset), lb_min(jset), &
     576              :                            tab, tkab, kab, rac, rbc, direction_Or)
     577              : 
     578              :                DEALLOCATE (tkab)
     579              :                ! Contraction step
     580              :                DO ir = 1, 3
     581              :                   CALL contraction(kab(:, :, ir), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
     582              :                                    cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), trans=trans)
     583              : 
     584              : !$                CALL omp_set_lock(locks(hash))
     585              :                   SELECT CASE (ir)
     586              :                   CASE (1)
     587              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kx_block, sgfa, sgfb, trans=trans)
     588              :                   CASE (2)
     589              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), ky_block, sgfa, sgfb, trans=trans)
     590              :                   CASE (3)
     591              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kz_block, sgfa, sgfb, trans=trans)
     592              :                   END SELECT
     593              : !$                CALL omp_unset_lock(locks(hash))
     594              :                END DO
     595              : 
     596              :             END DO
     597              :          END DO
     598              :       END DO !iterator
     599              :       DEALLOCATE (kab, qab)
     600              : !$OMP DO
     601              : !$    DO lock_num = 1, nlock
     602              : !$       call omp_destroy_lock(locks(lock_num))
     603              : !$    END DO
     604              : !$OMP END DO
     605              : 
     606              : !$OMP SINGLE
     607              : !$    DEALLOCATE (locks)
     608              : !$OMP END SINGLE NOWAIT
     609              : 
     610              : !$OMP END PARALLEL
     611              : 
     612              :       ! Release work storage
     613           76 :       DEALLOCATE (basis_set_list)
     614           76 :       CALL timestop(handle)
     615              : 
     616          228 :    END SUBROUTINE build_tr_matrix
     617              : 
     618              : ! **************************************************************************************************
     619              : !> \brief Commutator of the of the local part of the pseudopotential with r
     620              : !> \param matrix_rcore ...
     621              : !> \param qs_env ...
     622              : !> \param qs_kind_set ...
     623              : !> \param basis_type ...
     624              : !> \param sab_nl ...
     625              : !> \param rf ...
     626              : !> \author Edward Ditler, Tomas Zimmermann
     627              : ! **************************************************************************************************
     628           76 :    SUBROUTINE build_rcore_matrix(matrix_rcore, qs_env, qs_kind_set, basis_type, sab_nl, rf)
     629              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: matrix_rcore
     630              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     631              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     632              :       CHARACTER(LEN=*)                                   :: basis_type
     633              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     634              :          POINTER                                         :: sab_nl
     635              :       REAL(KIND=dp), DIMENSION(3), OPTIONAL              :: rf
     636              : 
     637              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_rcore_matrix'
     638              :       INTEGER, PARAMETER                                 :: nexp_max = 30
     639              : 
     640              :       INTEGER :: atom_a, atom_b, handle, i, iatom, icol, idir, ikind, inode, irow, iset, jatom, &
     641              :          jkind, jset, katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxlgto, maxlppl, maxnset, &
     642              :          maxsgf, mepos, n_local, ncoa, ncob, nder, nexp_lpot, nexp_ppl, nimages, nkind, nloc, &
     643              :          nseta, nsetb, nthread, sgfa, sgfb
     644           76 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     645              :       INTEGER, DIMENSION(1:10)                           :: nrloc
     646              :       INTEGER, DIMENSION(3)                              :: cellind
     647           76 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, &
     648           76 :                                                             nct_lpot, npgfa, npgfb, nsgfa, nsgfb
     649           76 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     650              :       INTEGER, DIMENSION(nexp_max)                       :: nct_ppl
     651              :       LOGICAL                                            :: do_symmetric, dokp, ecp_local, &
     652              :                                                             ecp_semi_local, found, lpotextended
     653              :       REAL(KIND=dp)                                      :: alpha, dab, dac, dbc, ppl_radius
     654           76 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: hab, qab
     655           76 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ppl_work, rhab, work
     656              :       REAL(KIND=dp), DIMENSION(1:10)                     :: aloc, bloc
     657              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, raf, rb, rbc, rbf
     658              :       REAL(KIND=dp), DIMENSION(4, nexp_max)              :: cval_ppl
     659           76 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_local, alpha_lpot, c_local, cexp_ppl, &
     660           76 :                                                             set_radius_a, set_radius_b
     661           76 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_lpot, hx_block, hy_block, hz_block, &
     662           76 :                                                             rpgfa, rpgfb, scon_a, scon_b, sphi_a, &
     663           76 :                                                             sphi_b, zeta, zetb
     664              :       REAL(KIND=dp), DIMENSION(nexp_max)                 :: alpha_ppl
     665           76 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     666              :       TYPE(cell_type), POINTER                           :: cell
     667              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     668           76 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     669              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     670              :       TYPE(neighbor_list_iterator_p_type), &
     671           76 :          DIMENSION(:), POINTER                           :: ap_iterator, nl_iterator
     672              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     673           76 :          POINTER                                         :: sac_ppl
     674           76 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     675              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
     676              : 
     677           76 :       CALL timeset(routineN, handle)
     678              : 
     679              :       CALL get_qs_env(qs_env=qs_env, &
     680              :                       atomic_kind_set=atomic_kind_set, &
     681              :                       qs_kind_set=qs_kind_set, &
     682              :                       particle_set=particle_set, &
     683              :                       sac_ppl=sac_ppl, &
     684           76 :                       cell=cell)
     685              : 
     686              :       ! check for symmetry
     687           76 :       CPASSERT(SIZE(sab_nl) > 0)
     688           76 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
     689              : 
     690           76 :       nkind = SIZE(qs_kind_set)
     691              : 
     692              :       ! prepare basis set
     693          380 :       ALLOCATE (basis_set_list(nkind))
     694           76 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
     695              : 
     696           76 :       nder = 0
     697           76 :       nimages = 1
     698              : 
     699              :       alpha_ppl = 0
     700              :       nct_ppl = 0
     701              :       cval_ppl = 0
     702              : 
     703           76 :       nkind = SIZE(atomic_kind_set)
     704              : 
     705           76 :       dokp = (nimages > 1)
     706              : 
     707           76 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     708              : 
     709           76 :       maxder = ncoset(nder)
     710              : 
     711              :       CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxlgto=maxlgto, &
     712              :                            maxsgf=maxsgf, maxnset=maxnset, maxlppl=maxlppl, &
     713           76 :                            basis_type=basis_type)
     714              : 
     715           76 :       maxl = MAX(maxlgto, maxlppl)
     716           76 :       CALL init_orbital_pointers(2*maxl + 2*nder + 2)
     717              : 
     718              :       !tz: maxco in maxco*ncoset(maxlgto+1) is an overkill,
     719              :       !    properly there should be maxpgf*ncoset(maxlgto+1), but maxpgf is difficult to get
     720           76 :       ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl, maxco*ncoset(maxlgto + 1))
     721           76 :       ldai = ncoset(2*maxlgto + 2)
     722              : 
     723          228 :       DO ikind = 1, nkind
     724          152 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
     725          228 :          IF (ASSOCIATED(basis_set_a)) THEN
     726          152 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     727              :          ELSE
     728            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     729              :          END IF
     730              :       END DO
     731              : 
     732              :       nthread = 1
     733           76 : !$    nthread = omp_get_max_threads()
     734              : 
     735           76 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
     736              : 
     737              :       ! iterator for basis/potential list
     738           76 :       CALL neighbor_list_iterator_create(ap_iterator, sac_ppl, search=.TRUE., nthread=nthread)
     739              : 
     740              : !$OMP PARALLEL &
     741              : !$OMP DEFAULT (NONE) &
     742              : !$OMP SHARED  (nl_iterator, ap_iterator, basis_set_list, &
     743              : !$OMP          atomic_kind_set, qs_kind_set, particle_set, &
     744              : !$OMP          sab_nl, sac_ppl, nthread, ncoset, nkind, &
     745              : !$OMP          atom_of_kind, ldsab,  maxnset, maxder, &
     746              : !$OMP          maxlgto, nder, maxco, dokp, cell) &
     747              : !$OMP SHARED (matrix_rcore, rf) &
     748              : !$OMP PRIVATE (ikind, jkind, inode, iatom, jatom, rab, basis_set_a, basis_set_b, atom_a) &
     749              : !$OMP PRIVATE (atom_b) &
     750              : !$OMP PRIVATE (nsetb) &
     751              : !$OMP PRIVATE (dab, irow, icol, hx_block, hy_block, hz_block, found, iset, ncoa) &
     752              : !$OMP PRIVATE (sgfa, jset, ncob, sgfb, work, hab, rhab, kkind, nseta) &
     753              : !$OMP PRIVATE (gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended) &
     754              : !$OMP PRIVATE (ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl) &
     755              : !$OMP PRIVATE (ecp_semi_local, nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc) &
     756              : !$OMP PRIVATE (mepos) &
     757              : !$OMP PRIVATE (katom, ppl_work, cellind, ecp_local) &
     758              : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, scon_a) &
     759              : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, sphi_b, zetb, scon_b) &
     760              : !$OMP PRIVATE (nloc, nrloc, aloc, bloc, n_local, a_local, c_local, ldai) &
     761           76 : !$OMP PRIVATE (ra, rb, qab, raf, rbf)
     762              : 
     763              :       mepos = 0
     764              : !$    mepos = omp_get_thread_num()
     765              : 
     766              :       ALLOCATE (hab(ldsab, ldsab), rhab(ldsab, ldsab, 3), work(ldsab, ldsab*(nder + 1), 3))
     767              :       ALLOCATE (qab(ldsab, ldsab))
     768              : 
     769              :       ldai = ncoset(2*maxlgto + 2)
     770              :       ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2 + 1)))
     771              : 
     772              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     773              : 
     774              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, inode=inode, &
     775              :                                 iatom=iatom, jatom=jatom, r=rab, cell=cellind)
     776              : 
     777              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     778              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     779              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     780              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     781              : 
     782              :          atom_a = atom_of_kind(iatom)
     783              :          atom_b = atom_of_kind(jatom)
     784              : 
     785              :          ! basis ikind
     786              :          first_sgfa => basis_set_a%first_sgf
     787              :          la_max => basis_set_a%lmax
     788              :          la_min => basis_set_a%lmin
     789              :          npgfa => basis_set_a%npgf
     790              :          nsgfa => basis_set_a%nsgf_set
     791              :          rpgfa => basis_set_a%pgf_radius
     792              :          set_radius_a => basis_set_a%set_radius
     793              :          sphi_a => basis_set_a%sphi
     794              :          zeta => basis_set_a%zet
     795              :          scon_a => basis_set_a%scon
     796              :          ! basis jkind
     797              :          first_sgfb => basis_set_b%first_sgf
     798              :          lb_max => basis_set_b%lmax
     799              :          lb_min => basis_set_b%lmin
     800              :          npgfb => basis_set_b%npgf
     801              :          nsgfb => basis_set_b%nsgf_set
     802              :          rpgfb => basis_set_b%pgf_radius
     803              :          set_radius_b => basis_set_b%set_radius
     804              :          sphi_b => basis_set_b%sphi
     805              :          zetb => basis_set_b%zet
     806              :          scon_b => basis_set_b%scon
     807              : 
     808              :          nseta = basis_set_a%nset
     809              :          nsetb = basis_set_b%nset
     810              : 
     811              :          ! *** Create matrix blocks for a new matrix block column ***
     812              :          irow = iatom
     813              :          icol = jatom
     814              : 
     815              :          NULLIFY (hx_block)
     816              :          NULLIFY (hy_block)
     817              :          NULLIFY (hz_block)
     818              :          CALL dbcsr_get_block_p(matrix=matrix_rcore(1)%matrix, &
     819              :                                 row=irow, col=icol, BLOCK=hx_block, found=found)
     820              :          CPASSERT(found)
     821              :          CALL dbcsr_get_block_p(matrix=matrix_rcore(2)%matrix, &
     822              :                                 row=irow, col=icol, BLOCK=hy_block, found=found)
     823              :          CPASSERT(found)
     824              :          CALL dbcsr_get_block_p(matrix=matrix_rcore(3)%matrix, &
     825              :                                 row=irow, col=icol, BLOCK=hz_block, found=found)
     826              :          CPASSERT(found)
     827              : 
     828              :          dab = NORM2(rab)
     829              :          ra = pbc(particle_set(iatom)%r(:), cell)
     830              :          rb(:) = ra(:) + rab(:)
     831              : 
     832              :          raf = pbc(rf, ra, cell)
     833              :          rbf = raf + rab
     834              : 
     835              :          ! loop over all kinds for pseudopotential atoms
     836              :          DO kkind = 1, nkind
     837              :             CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
     838              :                              sgp_potential=sgp_potential)
     839              :             IF (ASSOCIATED(gth_potential)) THEN
     840              :                CALL get_potential(potential=gth_potential, &
     841              :                                   alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
     842              :                                   lpot_present=lpotextended, ppl_radius=ppl_radius)
     843              :                nexp_ppl = 1
     844              :                alpha_ppl(1) = alpha
     845              :                nct_ppl(1) = SIZE(cexp_ppl)
     846              :                cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
     847              :                IF (lpotextended) THEN
     848              :                   CALL get_potential(potential=gth_potential, &
     849              :                                      nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, cval_lpot=cval_lpot)
     850              :                   CPASSERT(nexp_lpot < nexp_max)
     851              :                   nexp_ppl = nexp_lpot + 1
     852              :                   alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
     853              :                   nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
     854              :                   DO i = 1, nexp_lpot
     855              :                      cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
     856              :                   END DO
     857              :                END IF
     858              :             ELSE IF (ASSOCIATED(sgp_potential)) THEN
     859              :                CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
     860              :                                   ppl_radius=ppl_radius)
     861              :                IF (ecp_local) THEN
     862              :                   CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
     863              :                   IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
     864              :                   nexp_ppl = nloc
     865              :                   CPASSERT(nexp_ppl <= nexp_max)
     866              :                   nct_ppl(1:nloc) = nrloc(1:nloc) - 1
     867              :                   alpha_ppl(1:nloc) = bloc(1:nloc)
     868              :                   cval_ppl(1, 1:nloc) = aloc(1:nloc)
     869              :                ELSE
     870              :                   CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
     871              :                   nexp_ppl = n_local
     872              :                   CPASSERT(nexp_ppl <= nexp_max)
     873              :                   nct_ppl(1:n_local) = 1
     874              :                   alpha_ppl(1:n_local) = a_local(1:n_local)
     875              :                   cval_ppl(1, 1:n_local) = c_local(1:n_local)
     876              :                END IF
     877              :                IF (ecp_semi_local) THEN
     878              :                   CPABORT("VCD with semi-local ECPs not implemented")
     879              :                END IF
     880              :             ELSE
     881              :                CYCLE
     882              :             END IF
     883              : 
     884              :             CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
     885              : 
     886              :             DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
     887              : 
     888              :                CALL get_iterator_info(ap_iterator, mepos=mepos, jatom=katom, r=rac)
     889              :                dac = SQRT(SUM(rac*rac))
     890              :                rbc(:) = rac(:) - rab(:)
     891              :                dbc = SQRT(SUM(rbc*rbc))
     892              :                IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac) .OR. &
     893              :                    (MAXVAL(set_radius_b(:)) + ppl_radius < dbc)) THEN
     894              :                   CYCLE
     895              :                END IF
     896              :                DO iset = 1, nseta
     897              :                   IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
     898              :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     899              :                   !                  ncoa = npgfa(iset)*(ncoset(la_max(iset))-ncoset(la_min(iset)-1))
     900              :                   sgfa = first_sgfa(1, iset)
     901              :                   DO jset = 1, nsetb
     902              :                      IF (set_radius_b(jset) + ppl_radius < dbc) CYCLE
     903              :                      ncob = npgfb(jset)*ncoset(lb_max(jset))
     904              :                      !                     ncob = npgfb(jset)*(ncoset(lb_max(jset))-ncoset(lb_min(jset)-1))
     905              :                      sgfb = first_sgfb(1, jset)
     906              :                      IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     907              :                      ! *** Calculate the GTH pseudo potential forces ***
     908              :                      hab = 0
     909              :                      rhab = 0
     910              :                      ppl_work = 0
     911              :                      work = 0
     912              : 
     913              :                      CALL ppl_integral( &
     914              :                         la_max(iset) + 1, la_min(iset), npgfa(iset), &
     915              :                         rpgfa(:, iset), zeta(:, iset), &
     916              :                         lb_max(jset) + 1, lb_min(jset), npgfb(jset), &
     917              :                         rpgfb(:, jset), zetb(:, jset), &
     918              :                         nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
     919              :                         rab, dab, rac, dac, rbc, dbc, hab(:, :), ppl_work)
     920              : 
     921              :                      ! product with r
     922              :                      CALL ab_opr(la_max(iset), npgfa(iset), rpgfa(:, iset), 0, &
     923              :                                  lb_max(jset), npgfb(jset), rpgfb(:, jset), 0, &
     924              :                                  dab, hab(:, :), rhab(:, :, :), raf, rbf, &
     925              :                                  direction_Or=.FALSE.)
     926              : 
     927              :                      DO idir = 1, 3
     928              :                         CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     929              :                                    1.0_dp, rhab(1, 1, idir), SIZE(rhab, 1), &
     930              :                                    sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     931              :                                    0.0_dp, work(1, 1, idir), SIZE(work, 1))
     932              : !$OMP CRITICAL(h_block_critical)
     933              :                         SELECT CASE (idir)
     934              :                         CASE (1)
     935              :                            CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     936              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     937              :                                       work(1, 1, idir), SIZE(work, 1), &
     938              :                                       1.0_dp, hx_block(sgfa, sgfb), SIZE(hx_block, 1))
     939              :                         CASE (2)
     940              :                            CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     941              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     942              :                                       work(1, 1, idir), SIZE(work, 1), &
     943              :                                       1.0_dp, hy_block(sgfa, sgfb), SIZE(hy_block, 1))
     944              :                         CASE (3)
     945              :                            CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     946              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     947              :                                       work(1, 1, idir), SIZE(work, 1), &
     948              :                                       1.0_dp, hz_block(sgfa, sgfb), SIZE(hz_block, 1))
     949              :                         END SELECT
     950              : !$OMP END CRITICAL(h_block_critical)
     951              :                      END DO
     952              :                   END DO
     953              :                END DO
     954              :             END DO
     955              :          END DO
     956              :       END DO ! iterator
     957              : 
     958              :       DEALLOCATE (hab, rhab, work, ppl_work)
     959              : 
     960              : !$OMP END PARALLEL
     961              : 
     962           76 :       CALL neighbor_list_iterator_release(ap_iterator)
     963           76 :       CALL neighbor_list_iterator_release(nl_iterator)
     964              : 
     965           76 :       DEALLOCATE (atom_of_kind, basis_set_list)
     966              : 
     967           76 :       CALL timestop(handle)
     968              : 
     969          304 :    END SUBROUTINE build_rcore_matrix
     970              : 
     971              : ! **************************************************************************************************
     972              : !> \brief Commutator of the Hartree+XC potentials with r
     973              : !> \param matrix_rv ...
     974              : !> \param qs_env ...
     975              : !> \param rc ...
     976              : !> \author Edward Ditler, Tomas Zimmermann
     977              : ! **************************************************************************************************
     978           76 :    SUBROUTINE build_matrix_r_vhxc(matrix_rv, qs_env, rc)
     979              :       TYPE(dbcsr_p_type), DIMENSION(:, :), &
     980              :          INTENT(INOUT), POINTER                          :: matrix_rv
     981              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     982              :       REAL(KIND=dp), DIMENSION(3)                        :: rc
     983              : 
     984              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_matrix_r_vhxc'
     985              :       INTEGER, PARAMETER                                 :: nspins = 1
     986              : 
     987              :       INTEGER                                            :: handle, idir, ispin
     988              :       REAL(kind=dp)                                      :: edisp
     989           76 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     990           76 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_rvxc, matrix_rvxc_desymm
     991              :       TYPE(pw_env_type), POINTER                         :: pw_env
     992              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     993           76 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_rspace, v_tau_rspace
     994              :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     995              :       TYPE(qs_energy_type), POINTER                      :: energy
     996              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     997              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     998              :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     999              : 
    1000           76 :       CALL timeset(routineN, handle)
    1001              : 
    1002              :       CALL get_qs_env(qs_env, matrix_ks=matrix_ks, &
    1003              :                       ks_env=ks_env, &
    1004              :                       pw_env=pw_env, &
    1005              :                       input=input, &
    1006              :                       v_hartree_rspace=v_hartree_rspace, &
    1007           76 :                       energy=energy)
    1008              : 
    1009           76 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    1010              : 
    1011              :       ! matrix_rv(nspin, 3) was allocated as
    1012              :       ! CALL dbcsr_init_p(vcd_env%matrix_hr(ispin, i)%matrix)
    1013              :       ! CALL dbcsr_copy(vcd_env%matrix_hr(ispin, i)%matrix, vcd_env%matrix_nosym_temp(1)%matrix)
    1014              : 
    1015           76 :       NULLIFY (matrix_rvxc, matrix_rvxc_desymm)
    1016           76 :       CALL dbcsr_allocate_matrix_set(matrix_rvxc, nspins, 3)
    1017           76 :       CALL dbcsr_allocate_matrix_set(matrix_rvxc_desymm, nspins, 3)
    1018              : 
    1019          152 :       DO ispin = 1, nspins
    1020          380 :          DO idir = 1, 3
    1021          228 :             CALL dbcsr_init_p(matrix_rvxc(ispin, idir)%matrix)
    1022          228 :             CALL dbcsr_init_p(matrix_rvxc_desymm(ispin, idir)%matrix)
    1023              : 
    1024          228 :             CALL dbcsr_copy(matrix_rvxc_desymm(ispin, idir)%matrix, matrix_rv(1, 1)%matrix)
    1025          228 :             CALL dbcsr_set(matrix_rvxc_desymm(ispin, idir)%matrix, 0._dp)
    1026              : 
    1027          228 :             CALL dbcsr_copy(matrix_rvxc(ispin, idir)%matrix, matrix_ks(ispin)%matrix)
    1028          304 :             CALL dbcsr_set(matrix_rvxc(ispin, idir)%matrix, 0.0_dp)
    1029              :          END DO
    1030              :       END DO
    1031              : 
    1032           76 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
    1033           76 :       CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
    1034              : 
    1035           76 :       NULLIFY (v_rspace)
    1036           76 :       NULLIFY (v_tau_rspace)
    1037              :       CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
    1038              :                          vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=energy%exc, &
    1039              :                          edisp=edisp, dispersion_env=qs_env%dispersion_env, &
    1040           76 :                          just_energy=.FALSE.)
    1041              : 
    1042           76 :       IF (.NOT. ASSOCIATED(v_rspace)) THEN
    1043            0 :          ALLOCATE (v_rspace(nspins))
    1044            0 :          DO ispin = 1, nspins
    1045            0 :             CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
    1046            0 :             CALL pw_zero(v_rspace(ispin))
    1047              :          END DO
    1048              :       END IF
    1049              : 
    1050          152 :       DO ispin = 1, nspins
    1051           76 :          CALL pw_axpy(v_hartree_rspace, v_rspace(ispin), 1.0_dp/v_hartree_rspace%pw_grid%dvol)
    1052              :          CALL integrate_rv_rspace(v_rspace=v_rspace(ispin), hmat=matrix_rvxc(ispin, :), qs_env=qs_env, &
    1053           76 :                                   rc=rc)
    1054              : 
    1055          380 :          DO idir = 1, 3
    1056          228 :             CALL dbcsr_scale(matrix_rvxc(ispin, idir)%matrix, v_rspace(ispin)%pw_grid%dvol)
    1057          228 :             CALL dbcsr_desymmetrize(matrix_rvxc(ispin, idir)%matrix, matrix_rvxc_desymm(ispin, idir)%matrix)
    1058          304 :             CALL dbcsr_add(matrix_rv(ispin, idir)%matrix, matrix_rvxc_desymm(ispin, idir)%matrix, 1.0_dp, 1.0_dp)
    1059              :          END DO
    1060              :       END DO
    1061              : 
    1062              :       ! return pw grids
    1063          152 :       DO ispin = 1, nspins
    1064          152 :          CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
    1065              :       END DO
    1066           76 :       DEALLOCATE (v_rspace)
    1067              : 
    1068           76 :       CALL dbcsr_deallocate_matrix_set(matrix_rvxc)
    1069           76 :       CALL dbcsr_deallocate_matrix_set(matrix_rvxc_desymm)
    1070              : 
    1071           76 :       CALL timestop(handle)
    1072              : 
    1073           76 :    END SUBROUTINE build_matrix_r_vhxc
    1074              : 
    1075              : ! **************************************************************************************************
    1076              : !> \brief Calculates the integrals < mu | r * V | nu >
    1077              : !>        There is no direction_Or argument, because the potentials commute with r
    1078              : !>        This routine uses integrate_pgf_product directly. It could probably be rewritten to use
    1079              : !>          the new task_list interface.
    1080              : !> \param v_rspace ...
    1081              : !> \param hmat ...
    1082              : !> \param qs_env ...
    1083              : !> \param rc ...
    1084              : !> \author Edward Ditler
    1085              : ! **************************************************************************************************
    1086           76 :    SUBROUTINE integrate_rv_rspace(v_rspace, hmat, qs_env, rc)
    1087              : 
    1088              :       TYPE(pw_r3d_rs_type)                               :: v_rspace
    1089              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: hmat
    1090              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1091              :       REAL(KIND=dp), DIMENSION(3)                        :: rc
    1092              : 
    1093              :       CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rv_rspace'
    1094              : 
    1095              :       CHARACTER(len=default_string_length)               :: my_basis_type
    1096              :       INTEGER :: bcol, brow, handle, iatom, idir, igrid_level, ikind, ikind_old, ilevel, img, &
    1097              :          ipair, ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, &
    1098              :          jpgf, jpgf_new, jset, jset_new, jset_old, ldsab, maxco, maxlgto, maxpgf, maxset, &
    1099              :          maxsgf_set, na1, na2, natom, nb1, nb2, ncoa, ncoa_full, ncob, ncob_full, nkind, nseta, &
    1100              :          nsetb, nthread, sgfa, sgfb
    1101           76 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    1102           76 :                                                             npgfb, nsgfa, nsgfb
    1103           76 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    1104              :       LOGICAL :: atom_pair_changed, atom_pair_done, distributed_grids, found, has_threads, &
    1105              :          my_compute_tau, my_gapw, new_set_pair_coming
    1106              :       REAL(KIND=dp)                                      :: dab, eps_rho_rspace, f, prefactor, &
    1107              :                                                             radius, zetp
    1108              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rab_inv, rac, rb, rbc, rp
    1109           76 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
    1110           76 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, rpgfa, rpgfb, sphi_a, sphi_b, work, &
    1111           76 :                                                             zeta, zetb
    1112           76 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: habt, rhab, workt
    1113           76 :       TYPE(atom_pair_type), DIMENSION(:), POINTER        :: atom_pair_recv, atom_pair_send
    1114           76 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1115           76 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: h_block
    1116              :       TYPE(cell_type), POINTER                           :: cell
    1117              :       TYPE(dbcsr_distribution_type)                      :: dist
    1118              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1119              :       TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
    1120              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1121              :       TYPE(mp_comm_type)                                 :: group
    1122           76 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1123              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1124           76 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1125           76 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_v
    1126              :       TYPE(task_list_type), POINTER                      :: task_list, task_list_soft
    1127           76 :       TYPE(task_type), DIMENSION(:), POINTER             :: tasks
    1128              :       TYPE(virial_type), POINTER                         :: virial
    1129              : 
    1130           76 :       CALL timeset(routineN, handle)
    1131              : 
    1132           76 :       my_compute_tau = .FALSE.
    1133           76 :       my_gapw = .FALSE.
    1134           76 :       my_basis_type = "ORB"
    1135              : 
    1136              :       ! get the task lists
    1137              :       CALL get_qs_env(qs_env=qs_env, &
    1138              :                       task_list=task_list, &
    1139           76 :                       task_list_soft=task_list_soft)
    1140           76 :       CPASSERT(ASSOCIATED(task_list))
    1141              : 
    1142              :       ! the information on the grids is provided through pw_env
    1143              :       ! pw_env has to be the parent env for the potential grid (input)
    1144              :       ! there is an option to provide an external grid
    1145           76 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    1146           76 :       CPASSERT(ASSOCIATED(pw_env))
    1147              : 
    1148              :       ! get all the general information on the system we are working on
    1149              :       CALL get_qs_env(qs_env=qs_env, &
    1150              :                       atomic_kind_set=atomic_kind_set, &
    1151              :                       qs_kind_set=qs_kind_set, &
    1152              :                       cell=cell, &
    1153              :                       dft_control=dft_control, &
    1154              :                       particle_set=particle_set, &
    1155              :                       virial=virial, &
    1156           76 :                       natom=natom)
    1157              : 
    1158           76 :       rab = 0._dp
    1159           76 :       rac = 0._dp
    1160           76 :       rbc = 0._dp
    1161              : 
    1162              :       ! short cuts to task list variables
    1163           76 :       tasks => task_list%tasks
    1164           76 :       atom_pair_send => task_list%atom_pair_send
    1165           76 :       atom_pair_recv => task_list%atom_pair_recv
    1166              : 
    1167           76 :       CPASSERT(ASSOCIATED(pw_env))
    1168           76 :       CALL pw_env_get(pw_env, rs_grids=rs_v)
    1169              : 
    1170              :       ! get mpi group from rs_v
    1171           76 :       group = rs_v(1)%desc%group
    1172              : 
    1173              :       ! assign from pw_env
    1174           76 :       gridlevel_info => pw_env%gridlevel_info
    1175              : 
    1176              :       ! transform the potential on the rs_multigrids
    1177           76 :       CALL potential_pw2rs(rs_v, v_rspace, pw_env)
    1178              : 
    1179           76 :       nkind = SIZE(qs_kind_set)
    1180              : 
    1181              :       ! needs to be consistent with rho_rspace
    1182           76 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
    1183              : 
    1184              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
    1185              :                            maxco=maxco, &
    1186              :                            maxlgto=maxlgto, &
    1187              :                            maxsgf_set=maxsgf_set, &
    1188           76 :                            basis_type=my_basis_type)
    1189              : 
    1190           76 :       distributed_grids = .FALSE.
    1191          380 :       DO igrid_level = 1, gridlevel_info%ngrid_levels
    1192           76 :          IF (rs_v(igrid_level)%desc%distributed) THEN
    1193          304 :             distributed_grids = .TRUE.
    1194              :          END IF
    1195              :       END DO
    1196              : 
    1197              :       nthread = 1
    1198           76 : !$    nthread = omp_get_max_threads()
    1199              : 
    1200              :       ! get maximum numbers
    1201           76 :       maxset = 0
    1202           76 :       maxpgf = 0
    1203          228 :       DO ikind = 1, nkind
    1204              :          CALL get_qs_kind(qs_kind_set(ikind), &
    1205          152 :                           basis_set=orb_basis_set, basis_type=my_basis_type)
    1206              : 
    1207          152 :          IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
    1208              : 
    1209              :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1210          152 :                                 npgf=npgfa, nset=nseta)
    1211              : 
    1212          152 :          maxset = MAX(nseta, maxset)
    1213          380 :          maxpgf = MAX(MAXVAL(npgfa), maxpgf)
    1214              :       END DO
    1215              : 
    1216           76 :       ldsab = MAX(maxco, maxsgf_set, maxpgf*ncoset(maxlgto + 1))
    1217              : 
    1218              :       !   *** Allocate work storage ***
    1219           76 :       NULLIFY (habt, workt)
    1220           76 :       CALL reallocate(habt, 1, ldsab, 1, ldsab, 0, nthread)
    1221           76 :       CALL reallocate(workt, 1, ldsab, 1, maxsgf_set, 0, nthread)
    1222          380 :       ALLOCATE (rhab(ldsab, ldsab, 3))
    1223              : 
    1224          304 :       ALLOCATE (h_block(3))
    1225              : 
    1226           76 :       ithread = 0
    1227           76 : !$    ithread = omp_get_thread_num()
    1228           76 :       work => workt(:, :, ithread)
    1229           76 :       hab => habt(:, :, ithread)
    1230       124716 :       hab(:, :) = 0._dp
    1231              : 
    1232           76 :       iset_old = -1; jset_old = -1
    1233           76 :       ikind_old = -1; jkind_old = -1
    1234              : 
    1235              :       ! Here we loop over gridlevels first, finalising the matrix after each grid level is
    1236              :       ! completed.  On each grid level, we loop over atom pairs, which will only access
    1237              :       ! a single block of each matrix, so with OpenMP, each matrix block is only touched
    1238              :       ! by a single thread for each grid level
    1239          380 :       loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels
    1240         1216 :          DO idir = 1, 3
    1241          912 :             CALL dbcsr_work_create(hmat(idir)%matrix, work_mutable=.TRUE., n=nthread)
    1242          912 :             CALL dbcsr_get_info(hmat(idir)%matrix, distribution=dist)
    1243          912 :             CALL dbcsr_distribution_get(dist, has_threads=has_threads)
    1244          912 : !$          IF (.NOT. has_threads) &
    1245         1216 : !$             CPABORT("No thread distribution defined.")
    1246              :          END DO
    1247              : 
    1248          988 :          loop_pairs: DO ipair = 1, task_list%npairs(igrid_level)
    1249         4598 :          loop_tasks: DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
    1250         3610 :             ilevel = tasks(itask)%grid_level
    1251         3610 :             img = tasks(itask)%image
    1252         3610 :             iatom = tasks(itask)%iatom
    1253         3610 :             jatom = tasks(itask)%jatom
    1254         3610 :             iset = tasks(itask)%iset
    1255         3610 :             jset = tasks(itask)%jset
    1256         3610 :             ipgf = tasks(itask)%ipgf
    1257         3610 :             jpgf = tasks(itask)%jpgf
    1258         3610 :             CPASSERT(img == 1)
    1259              : 
    1260              :             ! At the start of a block of tasks, get atom data (and kind data, if needed)
    1261         3610 :             IF (itask .EQ. task_list%taskstart(ipair, igrid_level)) THEN
    1262              : 
    1263          684 :                ikind = particle_set(iatom)%atomic_kind%kind_number
    1264          684 :                jkind = particle_set(jatom)%atomic_kind%kind_number
    1265              : 
    1266          684 :                IF (iatom <= jatom) THEN
    1267          456 :                   brow = iatom
    1268          456 :                   bcol = jatom
    1269              :                ELSE
    1270          228 :                   brow = jatom
    1271          228 :                   bcol = iatom
    1272              :                END IF
    1273              : 
    1274          684 :                IF (ikind .NE. ikind_old) THEN
    1275              :                   CALL get_qs_kind(qs_kind_set(ikind), &
    1276           76 :                                    basis_set=orb_basis_set, basis_type=my_basis_type)
    1277              : 
    1278              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1279              :                                          first_sgf=first_sgfa, &
    1280              :                                          lmax=la_max, &
    1281              :                                          lmin=la_min, &
    1282              :                                          npgf=npgfa, &
    1283              :                                          nset=nseta, &
    1284              :                                          nsgf_set=nsgfa, &
    1285              :                                          pgf_radius=rpgfa, &
    1286              :                                          set_radius=set_radius_a, &
    1287              :                                          sphi=sphi_a, &
    1288           76 :                                          zet=zeta)
    1289              :                END IF
    1290              : 
    1291          684 :                IF (jkind .NE. jkind_old) THEN
    1292              :                   CALL get_qs_kind(qs_kind_set(jkind), &
    1293          456 :                                    basis_set=orb_basis_set, basis_type=my_basis_type)
    1294              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1295              :                                          first_sgf=first_sgfb, &
    1296              :                                          lmax=lb_max, &
    1297              :                                          lmin=lb_min, &
    1298              :                                          npgf=npgfb, &
    1299              :                                          nset=nsetb, &
    1300              :                                          nsgf_set=nsgfb, &
    1301              :                                          pgf_radius=rpgfb, &
    1302              :                                          set_radius=set_radius_b, &
    1303              :                                          sphi=sphi_b, &
    1304          456 :                                          zet=zetb)
    1305              : 
    1306              :                END IF
    1307              : 
    1308         2736 :                DO idir = 1, 3
    1309         2052 :                   NULLIFY (h_block(idir)%block)
    1310         2052 :                   CALL dbcsr_get_block_p(hmat(idir)%matrix, brow, bcol, h_block(idir)%block, found)
    1311         2736 :                   CPASSERT(found)
    1312              :                END DO
    1313              : 
    1314              :                ikind_old = ikind
    1315              :                jkind_old = jkind
    1316              : 
    1317              :                atom_pair_changed = .TRUE.
    1318              : 
    1319              :             ELSE
    1320              : 
    1321              :                atom_pair_changed = .FALSE.
    1322              : 
    1323              :             END IF
    1324              : 
    1325         3610 :             IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
    1326              :                ! We reuse the hab(:, :) array to put the new integrals in.
    1327              : 
    1328          684 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
    1329          684 :                ncoa_full = npgfa(iset)*ncoset(la_max(iset) + 1)
    1330          684 :                sgfa = first_sgfa(1, iset)
    1331          684 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
    1332          684 :                ncob_full = npgfb(jset)*ncoset(lb_max(jset) + 1)
    1333          684 :                sgfb = first_sgfb(1, jset)
    1334              : 
    1335          684 :                IF (iatom <= jatom) THEN
    1336       324216 :                   hab(1:ncoa_full, 1:ncob_full) = 0._dp
    1337              :                ELSE
    1338       106020 :                   hab(1:ncob_full, 1:ncoa_full) = 0._dp
    1339              :                END IF
    1340              : 
    1341              :                iset_old = iset
    1342              :                jset_old = jset
    1343              : 
    1344              :             END IF
    1345              : 
    1346        14440 :             rab = tasks(itask)%rab
    1347        14440 :             dab = norm2(rab)
    1348              :             ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
    1349         3610 :             ra = pbc(particle_set(iatom)%r(:), cell)
    1350        14440 :             rb(:) = ra(:) + rab(:)
    1351         3610 :             rac = pbc(rc, ra, cell)
    1352        14440 :             rbc = rac + rab
    1353              : 
    1354         3610 :             zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
    1355         3610 :             f = zetb(jpgf, jset)/zetp
    1356        14440 :             rp(:) = ra(:) + f*rab(:)
    1357              : 
    1358        14440 :             prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
    1359              :             radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
    1360              :                                               lb_min=lb_min(jset), lb_max=lb_max(jset), &
    1361              :                                               ra=ra, rb=rb, rp=rp, &
    1362              :                                               zetp=zetp, eps=eps_rho_rspace, &
    1363         3610 :                                               prefactor=prefactor, cutoff=1.0_dp)
    1364              : 
    1365         3610 :             na1 = (ipgf - 1)*ncoset(la_max(iset) + 1) + 1
    1366         3610 :             na2 = ipgf*ncoset(la_max(iset) + 1)
    1367         3610 :             nb1 = (jpgf - 1)*ncoset(lb_max(jset) + 1) + 1
    1368         3610 :             nb2 = jpgf*ncoset(lb_max(jset) + 1)
    1369              : 
    1370         3610 :             IF (iatom <= jatom) THEN
    1371              :                CALL integrate_pgf_product( &
    1372              :                   la_max(iset) + 1, zeta(ipgf, iset), la_min(iset), &
    1373              :                   lb_max(jset) + 1, zetb(jpgf, jset), lb_min(jset), &
    1374              :                   ra, rab, rs_v(igrid_level), &
    1375              :                   hab, o1=na1 - 1, o2=nb1 - 1, &
    1376              :                   radius=radius, &
    1377         2432 :                   calculate_forces=.FALSE.)
    1378              :             ELSE
    1379         4712 :                rab_inv = -rab
    1380              :                CALL integrate_pgf_product( &
    1381              :                   lb_max(jset) + 1, zetb(jpgf, jset), lb_min(jset), &
    1382              :                   la_max(iset) + 1, zeta(ipgf, iset), la_min(iset), &
    1383              :                   rb, rab_inv, rs_v(igrid_level), &
    1384              :                   hab, o1=nb1 - 1, o2=na1 - 1, &
    1385              :                   radius=radius, &
    1386         1178 :                   calculate_forces=.FALSE.)
    1387              :             END IF
    1388              : 
    1389         3610 :             new_set_pair_coming = .FALSE.
    1390         3610 :             atom_pair_done = .FALSE.
    1391         3610 :             IF (itask < task_list%taskstop(ipair, igrid_level)) THEN
    1392         2926 :                ilevel = tasks(itask + 1)%grid_level
    1393         2926 :                img = tasks(itask + 1)%image
    1394         2926 :                iatom = tasks(itask + 1)%iatom
    1395         2926 :                jatom = tasks(itask + 1)%jatom
    1396         2926 :                iset_new = tasks(itask + 1)%iset
    1397         2926 :                jset_new = tasks(itask + 1)%jset
    1398         2926 :                ipgf_new = tasks(itask + 1)%ipgf
    1399         2926 :                jpgf_new = tasks(itask + 1)%jpgf
    1400         2926 :                IF (iset_new .NE. iset .OR. jset_new .NE. jset) THEN
    1401            0 :                   new_set_pair_coming = .TRUE.
    1402              :                END IF
    1403              :             ELSE
    1404              :                ! do not forget the last block
    1405              :                new_set_pair_coming = .TRUE.
    1406         3610 :                atom_pair_done = .TRUE.
    1407              :             END IF
    1408              : 
    1409         4294 :             IF (new_set_pair_coming) THEN
    1410              :                ! Increase lx, ly, lz by one to account for the | r * b >
    1411          684 :                IF (iatom <= jatom) THEN
    1412              :                   ! direction_Or = .false. so that we use rac
    1413              :                   CALL ab_opr(la_max(iset), npgfa(iset), rpgfa(:, iset), 0, &
    1414              :                               lb_max(jset), npgfb(jset), rpgfb(:, jset), 0, &
    1415          456 :                               dab, hab(:, :), rhab(:, :, :), rac, rbc, direction_Or=.FALSE.)
    1416              : 
    1417              :                ELSE
    1418              :                   ! direction_Or = .true. so that we use rac
    1419              :                   CALL ab_opr(lb_max(jset), npgfb(jset), rpgfb(:, jset), 0, &
    1420              :                               la_max(iset), npgfa(iset), rpgfa(:, iset), 0, &
    1421          228 :                               dab, hab(:, :), rhab(:, :, :), rbc, rac, direction_Or=.TRUE.)
    1422              :                END IF
    1423              : 
    1424              :                ! contract the block into h if we're done with the current set pair
    1425         2736 :                DO idir = 1, 3
    1426         2736 :                   IF (iatom <= jatom) THEN
    1427       225720 :                      work = 0._dp
    1428       603630 :                      work(1:ncoa, 1:nsgfb(jset)) = MATMUL(rhab(1:ncoa, 1:ncob, idir), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
    1429              :                      h_block(idir)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
    1430              :                         h_block(idir)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
    1431       148428 :                         MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
    1432              :                   ELSE
    1433        72504 :                      work(1:ncob, 1:nsgfa(iset)) = MATMUL(rhab(1:ncob, 1:ncoa, idir), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
    1434              :                      h_block(idir)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
    1435              :                         h_block(idir)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
    1436        34200 :                         MATMUL(TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)), work(1:ncob, 1:nsgfa(iset)))
    1437              :                   END IF
    1438              :                END DO
    1439              :             END IF
    1440              : 
    1441              :          END DO loop_tasks
    1442              :          END DO loop_pairs
    1443              : 
    1444         1292 :          DO idir = 1, 3
    1445         1216 :             CALL dbcsr_finalize(hmat(idir)%matrix)
    1446              :          END DO
    1447              : 
    1448              :       END DO loop_gridlevels
    1449              : 
    1450              :       !   *** Release work storage ***
    1451           76 :       DEALLOCATE (habt, rhab, workt, h_block)
    1452              : 
    1453           76 :       CALL timestop(handle)
    1454              : 
    1455          152 :    END SUBROUTINE integrate_rv_rspace
    1456              : 
    1457              : ! **************************************************************************************************
    1458              : !> \brief Builds the overlap derivative wrt nuclear velocities
    1459              : !>         dS/dV = < mu | r | nu > * (nu - mu)
    1460              : !> \param qs_env ...
    1461              : !> \param matrix_dsdv ...
    1462              : !> \param deltaR ...
    1463              : !> \param rcc ...
    1464              : !> \author Edward Ditler
    1465              : ! **************************************************************************************************
    1466            6 :    SUBROUTINE build_dSdV_matrix(qs_env, matrix_dsdv, deltaR, rcc)
    1467              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1468              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: matrix_dsdv
    1469              :       REAL(KIND=dp), DIMENSION(:, :)                     :: deltaR
    1470              :       REAL(KIND=dp), DIMENSION(3)                        :: rcc
    1471              : 
    1472              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'build_dSdV_matrix'
    1473              : 
    1474              :       INTEGER                                            :: handle, i
    1475            6 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, my_matrix_dsdv, &
    1476            6 :                                                             my_matrix_dsdv2, my_matrix_mom
    1477              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1478            6 :          POINTER                                         :: sab_all
    1479            6 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1480              : 
    1481            6 :       CALL timeset(routineN, handle)
    1482              : 
    1483            6 :       NULLIFY (my_matrix_mom, my_matrix_dsdv, my_matrix_dsdv2, sab_all, qs_kind_set)
    1484              : 
    1485              :       CALL get_qs_env(qs_env=qs_env, &
    1486              :                       sab_all=sab_all, &
    1487              :                       qs_kind_set=qs_kind_set, &
    1488            6 :                       matrix_ks=matrix_ks)
    1489              : 
    1490              :       ! my_matrix_mom is needed because build_local_moment only works correctly with symmetric matrices
    1491            6 :       CALL dbcsr_allocate_matrix_set(my_matrix_dsdv, 3)
    1492            6 :       CALL dbcsr_allocate_matrix_set(my_matrix_dsdv2, 3)
    1493            6 :       CALL dbcsr_allocate_matrix_set(my_matrix_mom, 3)
    1494              : 
    1495           24 :       DO i = 1, 3
    1496           18 :          ALLOCATE (my_matrix_dsdv(i)%matrix)
    1497           18 :          ALLOCATE (my_matrix_dsdv2(i)%matrix)
    1498           18 :          ALLOCATE (my_matrix_mom(i)%matrix)
    1499              : 
    1500           18 :          CALL dbcsr_copy(my_matrix_dsdv(i)%matrix, matrix_dsdv(i)%matrix)
    1501           18 :          CALL dbcsr_copy(my_matrix_dsdv2(i)%matrix, matrix_dsdv(i)%matrix)
    1502           18 :          CALL dbcsr_copy(my_matrix_mom(i)%matrix, matrix_ks(1)%matrix)
    1503              : 
    1504           18 :          CALL dbcsr_set(my_matrix_dsdv2(i)%matrix, 0.0_dp)
    1505           18 :          CALL dbcsr_set(my_matrix_dsdv(i)%matrix, 0.0_dp)
    1506           18 :          CALL dbcsr_set(my_matrix_mom(i)%matrix, 0.0_dp)
    1507           24 :          CALL dbcsr_set(matrix_dsdv(i)%matrix, 0.0_dp)
    1508              :       END DO
    1509              : 
    1510            6 :       CALL build_dsdv_moments(qs_env, my_matrix_mom, 1, rcc)
    1511              : 
    1512           24 :       DO i = 1, 3
    1513           18 :          CALL dbcsr_desymmetrize(my_matrix_mom(i)%matrix, my_matrix_dsdv(i)%matrix)
    1514           24 :          CALL dbcsr_copy(my_matrix_dsdv2(i)%matrix, my_matrix_dsdv(i)%matrix)
    1515              :       END DO
    1516              : 
    1517              :       ! delta_nu^A <mu|r|nu>
    1518              :       CALL hr_mult_by_delta_3d(my_matrix_dsdv, qs_kind_set, "ORB", sab_all, &
    1519            6 :                                deltaR, direction_Or=.TRUE.)
    1520              :       ! -delta_mu^A <mu|r|nu>
    1521              :       CALL hr_mult_by_delta_3d(my_matrix_dsdv2, qs_kind_set, "ORB", sab_all, &
    1522            6 :                                deltaR, direction_Or=.FALSE.)
    1523           24 :       DO i = 1, 3
    1524           18 :          CALL dbcsr_copy(matrix_dsdv(i)%matrix, my_matrix_dsdv(i)%matrix)
    1525           24 :          CALL dbcsr_add(matrix_dsdv(i)%matrix, my_matrix_dsdv2(i)%matrix, 1.0_dp, -1.0_dp)
    1526              :       END DO
    1527              : 
    1528            6 :       CALL dbcsr_deallocate_matrix_set(my_matrix_dsdv)
    1529            6 :       CALL dbcsr_deallocate_matrix_set(my_matrix_dsdv2)
    1530            6 :       CALL dbcsr_deallocate_matrix_set(my_matrix_mom)
    1531              : 
    1532            6 :       CALL timestop(handle)
    1533              : 
    1534            6 :    END SUBROUTINE build_dSdV_matrix
    1535              : 
    1536              : ! **************************************************************************************************
    1537              : !> \brief Builds the  [Vnl, r] * r from either side
    1538              : !> \param matrix_rv ...
    1539              : !> \param qs_kind_set ...
    1540              : !> \param sab_all ...
    1541              : !> \param sap_ppnl ...
    1542              : !> \param eps_ppnl ...
    1543              : !> \param particle_set ...
    1544              : !> \param cell ...
    1545              : !> \param direction_Or ...
    1546              : !> \author Edward Ditler, Tomas Zimmermann
    1547              : ! **************************************************************************************************
    1548            4 :    SUBROUTINE build_com_rpnl_r(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, &
    1549              :                                particle_set, cell, direction_Or)
    1550              : 
    1551              :       TYPE(dbcsr_p_type), DIMENSION(:, :)                :: matrix_rv
    1552              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1553              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1554              :          POINTER                                         :: sab_all, sap_ppnl
    1555              :       REAL(KIND=dp), INTENT(IN)                          :: eps_ppnl
    1556              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
    1557              :          POINTER                                         :: particle_set
    1558              :       TYPE(cell_type), POINTER                           :: cell
    1559              :       LOGICAL, INTENT(IN)                                :: direction_Or
    1560              : 
    1561              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_com_rpnl_r'
    1562              : 
    1563              :       INTEGER                                            :: handle, i, iab, iac, iatom, ibc, icol, &
    1564              :                                                             ikind, irow, j, jatom, jkind, kac, &
    1565              :                                                             kbc, kkind, na, natom, nb, nkind, np, &
    1566              :                                                             slot
    1567              :       INTEGER, DIMENSION(3)                              :: cell_b
    1568              :       LOGICAL                                            :: found, ppnl_present
    1569              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    1570            4 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: achint, acint, bchint, bcint
    1571              :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
    1572           52 :       TYPE(block_p_type), DIMENSION(3, 3)                :: blocks_rvr
    1573            4 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
    1574              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1575            4 :       TYPE(sap_int_type), DIMENSION(:), POINTER          :: sap_int
    1576              : 
    1577              : !$    INTEGER(kind=omp_lock_kind), &
    1578            4 : !$       ALLOCATABLE, DIMENSION(:) :: locks
    1579              : !$    INTEGER                                            :: lock_num, hash
    1580              : !$    INTEGER, PARAMETER                                 :: nlock = 501
    1581              : 
    1582            4 :       ppnl_present = ASSOCIATED(sap_ppnl)
    1583            4 :       IF (.NOT. ppnl_present) RETURN
    1584              : 
    1585            4 :       CALL timeset(routineN, handle)
    1586            4 :       nkind = SIZE(qs_kind_set)
    1587            4 :       natom = SIZE(particle_set)
    1588              : 
    1589              :       ! sap_int needs to be shared as multiple threads need to access this
    1590            4 :       NULLIFY (sap_int)
    1591           28 :       ALLOCATE (sap_int(nkind*nkind))
    1592           20 :       DO i = 1, nkind*nkind
    1593           16 :          NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
    1594           20 :          sap_int(i)%nalist = 0
    1595              :       END DO
    1596              : 
    1597              :       ! We put zero as a reference point, because actually in the integrals we need two different ones:
    1598              :       !  < a | (r - R^\lambda_\beta) * [V, r_\alpha - R^\eta_\alpha] | b >
    1599              :       !  The first reference point can be added in a seperate step as the term will be
    1600              :       !        - R^\lambda_\beta * < a | [V, r_\alpha] | b >
    1601              :       !      = + R^\lambda_\beta * < a | [r_\alpha, V] | b >
    1602              :       !  The second reference point is not important, because it disappears in the commutator
    1603              :       CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=2, moment_mode=.TRUE., &
    1604            4 :                           particle_set=particle_set, cell=cell, refpoint=[0._dp, 0._dp, 0._dp])
    1605              : 
    1606              :       ! *** Set up a sorting index
    1607            4 :       CALL sap_sort(sap_int)
    1608              : 
    1609           20 :       ALLOCATE (basis_set(nkind))
    1610           12 :       DO ikind = 1, nkind
    1611            8 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
    1612           12 :          IF (ASSOCIATED(orb_basis_set)) THEN
    1613            8 :             basis_set(ikind)%gto_basis_set => orb_basis_set
    1614              :          ELSE
    1615            0 :             NULLIFY (basis_set(ikind)%gto_basis_set)
    1616              :          END IF
    1617              :       END DO
    1618              : 
    1619              :       ! *** All integrals needed have been calculated and stored in sap_int
    1620              :       ! *** We now calculate the commutator matrix elements
    1621              : 
    1622              : !$OMP PARALLEL &
    1623              : !$OMP DEFAULT (NONE) &
    1624              : !$OMP SHARED  (basis_set, matrix_rv, direction_Or, &
    1625              : !$OMP          sap_int, nkind, eps_ppnl, locks, sab_all) &
    1626              : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
    1627              : !$OMP          iab, irow, icol, blocks_rvr, lock_num, &
    1628              : !$OMP          found, iac, ibc, alist_ac, alist_bc, &
    1629              : !$OMP          na, np, nb, kkind, kac, kbc, i, j, &
    1630            4 : !$OMP          hash, natom, acint, bcint, achint, bchint)
    1631              : 
    1632              : !$OMP SINGLE
    1633              : !$    ALLOCATE (locks(nlock))
    1634              : !$OMP END SINGLE
    1635              : 
    1636              : !$OMP DO
    1637              : !$    DO lock_num = 1, nlock
    1638              : !$       call omp_init_lock(locks(lock_num))
    1639              : !$    END DO
    1640              : !$OMP END DO
    1641              : 
    1642              : !$OMP DO SCHEDULE(GUIDED)
    1643              : 
    1644              :       DO slot = 1, sab_all(1)%nl_size
    1645              : 
    1646              :          ikind = sab_all(1)%nlist_task(slot)%ikind
    1647              :          jkind = sab_all(1)%nlist_task(slot)%jkind
    1648              :          iatom = sab_all(1)%nlist_task(slot)%iatom
    1649              :          jatom = sab_all(1)%nlist_task(slot)%jatom
    1650              :          cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
    1651              :          rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
    1652              : 
    1653              :          IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
    1654              :          IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
    1655              :          iab = ikind + nkind*(jkind - 1)
    1656              : 
    1657              :          ! *** Create matrix blocks for a new matrix block column ***
    1658              :          irow = iatom
    1659              :          icol = jatom
    1660              :          DO i = 1, 3
    1661              :             DO j = 1, 3
    1662              :                ! t_alpha = MOD(i - 1, 3) + 1
    1663              :                ! t_beta = FLOOR(REAL(i - 1, dp)/3._dp) + 1
    1664              : 
    1665              :                CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, blocks_rvr(i, j)%block, found)
    1666              :                CPASSERT(found)
    1667              :             END DO
    1668              :          END DO
    1669              : 
    1670              :          ! loop over all kinds for projector atom
    1671              :          DO kkind = 1, nkind
    1672              :             iac = ikind + nkind*(kkind - 1)
    1673              :             ibc = jkind + nkind*(kkind - 1)
    1674              :             IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
    1675              :             IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
    1676              :             CALL get_alist(sap_int(iac), alist_ac, iatom)
    1677              :             CALL get_alist(sap_int(ibc), alist_bc, jatom)
    1678              :             IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
    1679              :             IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
    1680              :             DO kac = 1, alist_ac%nclist
    1681              :                DO kbc = 1, alist_bc%nclist
    1682              :                   IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
    1683              : 
    1684              :                   ! Some documention, considering the C1 O1 O2 molecule
    1685              :                   ! The integrals are <a|p>
    1686              :                   ! sap_int(1:2, 1:2) -> [(a=C, p=C), (a=C, p=O), (a=O, p=C), (a=O, p=0)]
    1687              :                   ! the (a=O, p=C) entry has an alist
    1688              :                   !  alist has two elements: O1, O2
    1689              :                   !     alist(O1) -> clist(C1)%acint are the integrals
    1690              :                   !     alist(O2) -> clist(C1)%acint are the integrals
    1691              : 
    1692              :                   IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
    1693              :                      IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
    1694              :                      acint => alist_ac%clist(kac)%acint
    1695              :                      bcint => alist_bc%clist(kbc)%acint
    1696              :                      achint => alist_ac%clist(kac)%achint
    1697              :                      bchint => alist_bc%clist(kbc)%achint
    1698              :                      na = SIZE(acint, 1)
    1699              :                      np = SIZE(acint, 2)
    1700              :                      nb = SIZE(bcint, 1)
    1701              : !$                   hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
    1702              : !$                   CALL omp_set_lock(locks(hash))
    1703              : 
    1704              :                      ! Template:
    1705              :                      ! blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
    1706              :                      !                                  MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xV
    1707              :                      ! 1, x, y, z, xx, xy, xz, yy, yz, zz
    1708              :                      IF (direction_Or) THEN
    1709              :                         ! (Vnl*r_alpha)*r_beta
    1710              : 
    1711              :                         ! This part is symmetric because [r_beta, r_alpha] = 0
    1712              :                         ! matrix_rvr(x, gamma) += <mu|p>h_ij * <p|x*gamma|nu>
    1713              :                         blocks_rvr(1, 1)%block(1:na, 1:nb) = blocks_rvr(1, 1)%block(1:na, 1:nb) + &
    1714              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xx)))
    1715              :                         blocks_rvr(1, 2)%block(1:na, 1:nb) = blocks_rvr(1, 2)%block(1:na, 1:nb) + &
    1716              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
    1717              :                         blocks_rvr(1, 3)%block(1:na, 1:nb) = blocks_rvr(1, 3)%block(1:na, 1:nb) + &
    1718              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
    1719              : 
    1720              :                         ! matrix_rvr(y, gamma) += <mu|p>h_ij * <p|y*gamma|nu>
    1721              :                         blocks_rvr(2, 1)%block(1:na, 1:nb) = blocks_rvr(2, 1)%block(1:na, 1:nb) + &
    1722              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
    1723              :                         blocks_rvr(2, 2)%block(1:na, 1:nb) = blocks_rvr(2, 2)%block(1:na, 1:nb) + &
    1724              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yy)))
    1725              :                         blocks_rvr(2, 3)%block(1:na, 1:nb) = blocks_rvr(2, 3)%block(1:na, 1:nb) + &
    1726              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
    1727              : 
    1728              :                         ! matrix_rvr(z, gamma) += <mu|p>h_ij * <p|z*gamma|nu>
    1729              :                         blocks_rvr(3, 1)%block(1:na, 1:nb) = blocks_rvr(3, 1)%block(1:na, 1:nb) + &
    1730              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
    1731              :                         blocks_rvr(3, 2)%block(1:na, 1:nb) = blocks_rvr(3, 2)%block(1:na, 1:nb) + &
    1732              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
    1733              :                         blocks_rvr(3, 3)%block(1:na, 1:nb) = blocks_rvr(3, 3)%block(1:na, 1:nb) + &
    1734              :                                                              MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_zz)))
    1735              : 
    1736              :                         ! -(r_alpha*Vnl)*r_beta
    1737              :                         ! matrix_rvr(x, gamma) += <mu|gamma|p>h_ij * <p|x|nu>
    1738              :                         blocks_rvr(1, 1)%block(1:na, 1:nb) = blocks_rvr(1, 1)%block(1:na, 1:nb) - &
    1739              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1740              :                         blocks_rvr(1, 2)%block(1:na, 1:nb) = blocks_rvr(1, 2)%block(1:na, 1:nb) - &
    1741              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1742              :                         blocks_rvr(1, 3)%block(1:na, 1:nb) = blocks_rvr(1, 3)%block(1:na, 1:nb) - &
    1743              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1744              : 
    1745              :                         ! matrix_rvr(y, gamma) += <mu|gamma|p>h_ij * <p|y|nu>
    1746              :                         blocks_rvr(2, 1)%block(1:na, 1:nb) = blocks_rvr(2, 1)%block(1:na, 1:nb) - &
    1747              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1748              :                         blocks_rvr(2, 2)%block(1:na, 1:nb) = blocks_rvr(2, 2)%block(1:na, 1:nb) - &
    1749              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1750              :                         blocks_rvr(2, 3)%block(1:na, 1:nb) = blocks_rvr(2, 3)%block(1:na, 1:nb) - &
    1751              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1752              : 
    1753              :                         ! matrix_rvr(z, gamma) += <mu|gamma|p>h_ij * <p|z|nu>
    1754              :                         blocks_rvr(3, 1)%block(1:na, 1:nb) = blocks_rvr(3, 1)%block(1:na, 1:nb) - &
    1755              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1756              :                         blocks_rvr(3, 2)%block(1:na, 1:nb) = blocks_rvr(3, 2)%block(1:na, 1:nb) - &
    1757              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1758              :                         blocks_rvr(3, 3)%block(1:na, 1:nb) = blocks_rvr(3, 3)%block(1:na, 1:nb) - &
    1759              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1760              : 
    1761              :                         ! This means that we have
    1762              :                         ! blocks_rvr(1:3, 1:3) = blocks_rvr(alpha, beta)
    1763              :                      ELSE
    1764              :                         ! r_beta*(Vnl*r_alpha)
    1765              :                         blocks_rvr(1, 1)%block(1:na, 1:nb) = blocks_rvr(1, 1)%block(1:na, 1:nb) + &
    1766              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1767              :                         blocks_rvr(1, 2)%block(1:na, 1:nb) = blocks_rvr(1, 2)%block(1:na, 1:nb) + &
    1768              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1769              :                         blocks_rvr(1, 3)%block(1:na, 1:nb) = blocks_rvr(1, 3)%block(1:na, 1:nb) + &
    1770              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    1771              : 
    1772              :                         blocks_rvr(2, 1)%block(1:na, 1:nb) = blocks_rvr(2, 1)%block(1:na, 1:nb) + &
    1773              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1774              :                         blocks_rvr(2, 2)%block(1:na, 1:nb) = blocks_rvr(2, 2)%block(1:na, 1:nb) + &
    1775              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1776              :                         blocks_rvr(2, 3)%block(1:na, 1:nb) = blocks_rvr(2, 3)%block(1:na, 1:nb) + &
    1777              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    1778              : 
    1779              :                         blocks_rvr(3, 1)%block(1:na, 1:nb) = blocks_rvr(3, 1)%block(1:na, 1:nb) + &
    1780              :                                                              MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1781              :                         blocks_rvr(3, 2)%block(1:na, 1:nb) = blocks_rvr(3, 2)%block(1:na, 1:nb) + &
    1782              :                                                              MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1783              :                         blocks_rvr(3, 3)%block(1:na, 1:nb) = blocks_rvr(3, 3)%block(1:na, 1:nb) + &
    1784              :                                                              MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    1785              : 
    1786              :                         ! -r_beta*(r_alpha*Vnl)
    1787              :                         blocks_rvr(1, 1)%block(1:na, 1:nb) = blocks_rvr(1, 1)%block(1:na, 1:nb) - &
    1788              :                                                              MATMUL(achint(1:na, 1:np, bi_xx), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1789              :                         blocks_rvr(1, 2)%block(1:na, 1:nb) = blocks_rvr(1, 2)%block(1:na, 1:nb) - &
    1790              :                                                              MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1791              :                         blocks_rvr(1, 3)%block(1:na, 1:nb) = blocks_rvr(1, 3)%block(1:na, 1:nb) - &
    1792              :                                                              MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1793              : 
    1794              :                         blocks_rvr(2, 1)%block(1:na, 1:nb) = blocks_rvr(2, 1)%block(1:na, 1:nb) - &
    1795              :                                                              MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1796              :                         blocks_rvr(2, 2)%block(1:na, 1:nb) = blocks_rvr(2, 2)%block(1:na, 1:nb) - &
    1797              :                                                              MATMUL(achint(1:na, 1:np, bi_yy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1798              :                         blocks_rvr(2, 3)%block(1:na, 1:nb) = blocks_rvr(2, 3)%block(1:na, 1:nb) - &
    1799              :                                                              MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1800              : 
    1801              :                         blocks_rvr(3, 1)%block(1:na, 1:nb) = blocks_rvr(3, 1)%block(1:na, 1:nb) - &
    1802              :                                                              MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1803              :                         blocks_rvr(3, 2)%block(1:na, 1:nb) = blocks_rvr(3, 2)%block(1:na, 1:nb) - &
    1804              :                                                              MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1805              :                         blocks_rvr(3, 3)%block(1:na, 1:nb) = blocks_rvr(3, 3)%block(1:na, 1:nb) - &
    1806              :                                                              MATMUL(achint(1:na, 1:np, bi_zz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    1807              :                      END IF
    1808              : 
    1809              : !$                   CALL omp_unset_lock(locks(hash))
    1810              :                      EXIT ! We have found a match and there can be only one single match
    1811              :                   END IF
    1812              :                END DO
    1813              :             END DO
    1814              :          END DO
    1815              :          DO i = 1, 3
    1816              :             NULLIFY (blocks_rvr(i, 1)%block)
    1817              :             NULLIFY (blocks_rvr(i, 2)%block)
    1818              :             NULLIFY (blocks_rvr(i, 3)%block)
    1819              :          END DO
    1820              :       END DO
    1821              : 
    1822              : !$OMP DO
    1823              : !$    DO lock_num = 1, nlock
    1824              : !$       call omp_destroy_lock(locks(lock_num))
    1825              : !$    END DO
    1826              : !$OMP END DO
    1827              : 
    1828              : !$OMP SINGLE
    1829              : !$    DEALLOCATE (locks)
    1830              : !$OMP END SINGLE NOWAIT
    1831              : 
    1832              : !$OMP END PARALLEL
    1833              : 
    1834            4 :       CALL release_sap_int(sap_int)
    1835              : 
    1836            4 :       DEALLOCATE (basis_set)
    1837              : 
    1838            4 :       CALL timestop(handle)
    1839              : 
    1840           12 :    END SUBROUTINE build_com_rpnl_r
    1841              : 
    1842              : ! **************************************************************************************************
    1843              : !> \brief Calculate the double commutator [[Vnl, r], r]
    1844              : !> \param matrix_rv ...
    1845              : !> \param qs_kind_set ...
    1846              : !> \param sab_orb ...
    1847              : !> \param sap_ppnl ...
    1848              : !> \param eps_ppnl ...
    1849              : !> \param particle_set ...
    1850              : !> \param pseudoatom Only consider pseudopotentials on atom lambda
    1851              : !> \author Edward Ditler
    1852              : ! **************************************************************************************************
    1853            6 :    SUBROUTINE build_dcom_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, pseudoatom)
    1854              : 
    1855              :       TYPE(dbcsr_p_type), DIMENSION(:, :)                :: matrix_rv
    1856              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1857              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1858              :          POINTER                                         :: sab_orb, sap_ppnl
    1859              :       REAL(KIND=dp), INTENT(IN)                          :: eps_ppnl
    1860              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
    1861              :          POINTER                                         :: particle_set
    1862              :       INTEGER, INTENT(IN)                                :: pseudoatom
    1863              : 
    1864              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_dcom_rpnl'
    1865              : 
    1866              :       INTEGER                                            :: handle, i, iab, iac, iatom, ibc, icol, &
    1867              :                                                             ikind, irow, j, jatom, jkind, kac, &
    1868              :                                                             kbc, kkind, na, natom, nb, nkind, np, &
    1869              :                                                             slot
    1870              :       INTEGER, DIMENSION(3)                              :: cell_b
    1871              :       LOGICAL                                            :: found, ppnl_present
    1872              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    1873            6 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: achint, acint, bchint, bcint
    1874              :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
    1875           78 :       TYPE(block_p_type), DIMENSION(3, 3)                :: blocks_rv
    1876            6 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
    1877              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1878            6 :       TYPE(sap_int_type), DIMENSION(:), POINTER          :: sap_int
    1879              : 
    1880              : !$    INTEGER(kind=omp_lock_kind), &
    1881            6 : !$       ALLOCATABLE, DIMENSION(:) :: locks
    1882              : !$    INTEGER                                            :: lock_num, hash
    1883              : !$    INTEGER, PARAMETER                                 :: nlock = 501
    1884              : 
    1885            6 :       ppnl_present = ASSOCIATED(sap_ppnl)
    1886            6 :       IF (.NOT. ppnl_present) RETURN
    1887              : 
    1888            6 :       CALL timeset(routineN, handle)
    1889            6 :       nkind = SIZE(qs_kind_set)
    1890            6 :       natom = SIZE(particle_set)
    1891              : 
    1892              :       ! sap_int needs to be shared as multiple threads need to access this
    1893            6 :       NULLIFY (sap_int)
    1894           42 :       ALLOCATE (sap_int(nkind*nkind))
    1895           30 :       DO i = 1, nkind*nkind
    1896           24 :          NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
    1897           30 :          sap_int(i)%nalist = 0
    1898              :       END DO
    1899              : 
    1900              :       ! "nder" in moment_mode is "order"
    1901              :       CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=2, moment_mode=.TRUE., &
    1902            6 :                           particle_set=particle_set, pseudoatom=pseudoatom)
    1903              : 
    1904              :       ! *** Set up a sorting index
    1905            6 :       CALL sap_sort(sap_int)
    1906              : 
    1907           30 :       ALLOCATE (basis_set(nkind))
    1908           18 :       DO ikind = 1, nkind
    1909           12 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
    1910           18 :          IF (ASSOCIATED(orb_basis_set)) THEN
    1911           12 :             basis_set(ikind)%gto_basis_set => orb_basis_set
    1912              :          ELSE
    1913            0 :             NULLIFY (basis_set(ikind)%gto_basis_set)
    1914              :          END IF
    1915              :       END DO
    1916              : 
    1917              :       ! *** All integrals needed have been calculated and stored in sap_int
    1918              :       ! *** We now calculate the commutator matrix elements
    1919              : 
    1920              : !$OMP PARALLEL &
    1921              : !$OMP DEFAULT (NONE) &
    1922              : !$OMP SHARED  (basis_set, matrix_rv, &
    1923              : !$OMP          sap_int, nkind, eps_ppnl, locks, sab_orb) &
    1924              : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
    1925              : !$OMP          iab, irow, icol, blocks_rv, lock_num, &
    1926              : !$OMP          found, iac, ibc, alist_ac, alist_bc, &
    1927              : !$OMP          na, np, nb, kkind, kac, kbc, i, j, &
    1928            6 : !$OMP          hash, natom, acint, bcint, achint, bchint)
    1929              : 
    1930              : !$OMP SINGLE
    1931              : !$    ALLOCATE (locks(nlock))
    1932              : !$OMP END SINGLE
    1933              : 
    1934              : !$OMP DO
    1935              : !$    DO lock_num = 1, nlock
    1936              : !$       call omp_init_lock(locks(lock_num))
    1937              : !$    END DO
    1938              : !$OMP END DO
    1939              : 
    1940              : !$OMP DO SCHEDULE(GUIDED)
    1941              : 
    1942              :       DO slot = 1, sab_orb(1)%nl_size
    1943              : 
    1944              :          ikind = sab_orb(1)%nlist_task(slot)%ikind
    1945              :          jkind = sab_orb(1)%nlist_task(slot)%jkind
    1946              :          iatom = sab_orb(1)%nlist_task(slot)%iatom
    1947              :          jatom = sab_orb(1)%nlist_task(slot)%jatom
    1948              :          cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
    1949              :          rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
    1950              : 
    1951              :          IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
    1952              :          IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
    1953              :          iab = ikind + nkind*(jkind - 1)
    1954              : 
    1955              :          ! *** Create matrix blocks for a new matrix block column ***
    1956              :          IF (iatom <= jatom) THEN
    1957              :             irow = iatom
    1958              :             icol = jatom
    1959              :          ELSE
    1960              :             irow = jatom
    1961              :             icol = iatom
    1962              :          END IF
    1963              : 
    1964              :          DO i = 1, 3
    1965              :             DO j = 1, 3
    1966              :                CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, blocks_rv(i, j)%block, found)
    1967              :                blocks_rv(i, j)%block = 0._dp
    1968              :                CPASSERT(found)
    1969              :             END DO
    1970              :          END DO
    1971              : 
    1972              :          ! loop over all kinds for projector atom
    1973              :          DO kkind = 1, nkind
    1974              :             iac = ikind + nkind*(kkind - 1)
    1975              :             ibc = jkind + nkind*(kkind - 1)
    1976              :             IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
    1977              :             IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
    1978              :             CALL get_alist(sap_int(iac), alist_ac, iatom)
    1979              :             CALL get_alist(sap_int(ibc), alist_bc, jatom)
    1980              :             IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
    1981              :             IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
    1982              :             DO kac = 1, alist_ac%nclist
    1983              :                DO kbc = 1, alist_bc%nclist
    1984              :                   IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
    1985              :                   IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
    1986              :                      IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
    1987              :                      acint => alist_ac%clist(kac)%acint
    1988              :                      bcint => alist_bc%clist(kbc)%acint
    1989              :                      achint => alist_ac%clist(kac)%achint
    1990              :                      bchint => alist_bc%clist(kbc)%achint
    1991              :                      na = SIZE(acint, 1)
    1992              :                      np = SIZE(acint, 2)
    1993              :                      nb = SIZE(bcint, 1)
    1994              : !$                   hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
    1995              : !$                   CALL omp_set_lock(locks(hash))
    1996              :                      ! Template:
    1997              :                      ! blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
    1998              :                      !                                  MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xV
    1999              :                      IF (iatom <= jatom) THEN
    2000              :                         ! r_alpha*Vnl*r_beta
    2001              :                         blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) + &
    2002              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2003              : 
    2004              :                         blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) + &
    2005              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2006              : 
    2007              :                         blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) + &
    2008              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2009              : 
    2010              :                         blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) + &
    2011              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2012              : 
    2013              :                         blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) + &
    2014              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2015              : 
    2016              :                         blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) + &
    2017              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2018              : 
    2019              :                         blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) + &
    2020              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2021              : 
    2022              :                         blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) + &
    2023              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2024              : 
    2025              :                         blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) + &
    2026              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2027              : 
    2028              :                         ! -r_alpha*r_beta*Vnl
    2029              :                         blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) - &
    2030              :                                                             MATMUL(achint(1:na, 1:np, bi_xx), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2031              : 
    2032              :                         blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) - &
    2033              :                                                             MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2034              : 
    2035              :                         blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) - &
    2036              :                                                             MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2037              : 
    2038              :                         blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) - &
    2039              :                                                             MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2040              : 
    2041              :                         blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) - &
    2042              :                                                             MATMUL(achint(1:na, 1:np, bi_yy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2043              : 
    2044              :                         blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) - &
    2045              :                                                             MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2046              : 
    2047              :                         blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) - &
    2048              :                                                             MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2049              : 
    2050              :                         blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) - &
    2051              :                                                             MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2052              : 
    2053              :                         blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) - &
    2054              :                                                             MATMUL(achint(1:na, 1:np, bi_zz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2055              : 
    2056              :                         ! -Vnl*r_beta*r_alpha
    2057              :                         blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) - &
    2058              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xx)))
    2059              : 
    2060              :                         blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) - &
    2061              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
    2062              : 
    2063              :                         blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) - &
    2064              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
    2065              : 
    2066              :                         blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) - &
    2067              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
    2068              : 
    2069              :                         blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) - &
    2070              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yy)))
    2071              : 
    2072              :                         blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) - &
    2073              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
    2074              : 
    2075              :                         blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) - &
    2076              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
    2077              : 
    2078              :                         blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) - &
    2079              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
    2080              : 
    2081              :                         blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) - &
    2082              :                                                             MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_zz)))
    2083              : 
    2084              :                         ! +r_beta*Vnl*r_alpha
    2085              :                         blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) + &
    2086              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2087              : 
    2088              :                         blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) + &
    2089              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2090              : 
    2091              :                         blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) + &
    2092              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2093              : 
    2094              :                         blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) + &
    2095              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2096              : 
    2097              :                         blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) + &
    2098              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2099              : 
    2100              :                         blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) + &
    2101              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2102              : 
    2103              :                         blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) + &
    2104              :                                                             MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2105              : 
    2106              :                         blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) + &
    2107              :                                                             MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2108              : 
    2109              :                         blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) + &
    2110              :                                                             MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2111              :                      ELSE
    2112              :                         ! r_alpha*Vnl*r_beta
    2113              :                         blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) + &
    2114              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2115              : 
    2116              :                         blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) + &
    2117              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2118              : 
    2119              :                         blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) + &
    2120              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2121              : 
    2122              :                         blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) + &
    2123              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2124              : 
    2125              :                         blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) + &
    2126              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2127              : 
    2128              :                         blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) + &
    2129              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2130              : 
    2131              :                         blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) + &
    2132              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2133              : 
    2134              :                         blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) + &
    2135              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2136              : 
    2137              :                         blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) + &
    2138              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2139              : 
    2140              :                         ! -r_alpha*r_beta*Vnl
    2141              :                         blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) - &
    2142              :                                                             MATMUL(bchint(1:nb, 1:np, bi_xx), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2143              : 
    2144              :                         blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) - &
    2145              :                                                             MATMUL(bchint(1:nb, 1:np, bi_xy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2146              : 
    2147              :                         blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) - &
    2148              :                                                             MATMUL(bchint(1:nb, 1:np, bi_xz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2149              : 
    2150              :                         blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) - &
    2151              :                                                             MATMUL(bchint(1:nb, 1:np, bi_xy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2152              : 
    2153              :                         blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) - &
    2154              :                                                             MATMUL(bchint(1:nb, 1:np, bi_yy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2155              : 
    2156              :                         blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) - &
    2157              :                                                             MATMUL(bchint(1:nb, 1:np, bi_yz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2158              : 
    2159              :                         blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) - &
    2160              :                                                             MATMUL(bchint(1:nb, 1:np, bi_xz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2161              : 
    2162              :                         blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) - &
    2163              :                                                             MATMUL(bchint(1:nb, 1:np, bi_yz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2164              : 
    2165              :                         blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) - &
    2166              :                                                             MATMUL(bchint(1:nb, 1:np, bi_zz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
    2167              : 
    2168              :                         ! -Vnl*r_beta*r_alpha
    2169              :                         blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) - &
    2170              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xx)))
    2171              : 
    2172              :                         blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) - &
    2173              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xy)))
    2174              : 
    2175              :                         blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) - &
    2176              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xz)))
    2177              : 
    2178              :                         blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) - &
    2179              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xy)))
    2180              : 
    2181              :                         blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) - &
    2182              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yy)))
    2183              : 
    2184              :                         blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) - &
    2185              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yz)))
    2186              : 
    2187              :                         blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) - &
    2188              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xz)))
    2189              : 
    2190              :                         blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) - &
    2191              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yz)))
    2192              : 
    2193              :                         blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) - &
    2194              :                                                             MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_zz)))
    2195              : 
    2196              :                         ! +r_beta*Vnl*r_alpha
    2197              :                         blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) + &
    2198              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2199              : 
    2200              :                         blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) + &
    2201              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2202              : 
    2203              :                         blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) + &
    2204              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_x)))
    2205              : 
    2206              :                         blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) + &
    2207              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2208              : 
    2209              :                         blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) + &
    2210              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2211              : 
    2212              :                         blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) + &
    2213              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_y)))
    2214              : 
    2215              :                         blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) + &
    2216              :                                                             MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2217              : 
    2218              :                         blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) + &
    2219              :                                                             MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2220              : 
    2221              :                         blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) + &
    2222              :                                                             MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_z)))
    2223              : 
    2224              :                      END IF
    2225              : !$                   CALL omp_unset_lock(locks(hash))
    2226              :                      EXIT ! We have found a match and there can be only one single match
    2227              :                   END IF
    2228              :                END DO
    2229              :             END DO
    2230              :          END DO
    2231              :          DO i = 1, 3
    2232              :             NULLIFY (blocks_rv(i, 1)%block)
    2233              :             NULLIFY (blocks_rv(i, 2)%block)
    2234              :             NULLIFY (blocks_rv(i, 3)%block)
    2235              :          END DO
    2236              :       END DO
    2237              : 
    2238              : !$OMP DO
    2239              : !$    DO lock_num = 1, nlock
    2240              : !$       call omp_destroy_lock(locks(lock_num))
    2241              : !$    END DO
    2242              : !$OMP END DO
    2243              : 
    2244              : !$OMP SINGLE
    2245              : !$    DEALLOCATE (locks)
    2246              : !$OMP END SINGLE NOWAIT
    2247              : 
    2248              : !$OMP END PARALLEL
    2249              : 
    2250            6 :       CALL release_sap_int(sap_int)
    2251              : 
    2252            6 :       DEALLOCATE (basis_set)
    2253              : 
    2254            6 :       CALL timestop(handle)
    2255              : 
    2256           18 :    END SUBROUTINE build_dcom_rpnl
    2257              : 
    2258              : ! **************************************************************************************************
    2259              : !> \brief dV_nl/dV. Adapted from build_com_rpnl.
    2260              : !> \param matrix_rv ...
    2261              : !> \param qs_kind_set ...
    2262              : !> \param sab_all ...
    2263              : !> \param sap_ppnl ...
    2264              : !> \param eps_ppnl ...
    2265              : !> \param particle_set ...
    2266              : !> \param pseudoatom ...
    2267              : !> \author Edward Ditler
    2268              : ! **************************************************************************************************
    2269            6 :    SUBROUTINE build_drpnl_matrix(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, pseudoatom)
    2270              : 
    2271              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: matrix_rv
    2272              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2273              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2274              :          POINTER                                         :: sab_all, sap_ppnl
    2275              :       REAL(KIND=dp), INTENT(IN)                          :: eps_ppnl
    2276              :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
    2277              :          POINTER                                         :: particle_set
    2278              :       INTEGER                                            :: pseudoatom
    2279              : 
    2280              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_drpnl_matrix'
    2281              : 
    2282              :       INTEGER                                            :: handle, i, iab, iac, iatom, ibc, icol, &
    2283              :                                                             ikind, irow, jatom, jkind, kac, kbc, &
    2284              :                                                             kkind, na, natom, nb, nkind, np, slot
    2285              :       INTEGER, DIMENSION(3)                              :: cell_b
    2286              :       LOGICAL                                            :: found, ppnl_present
    2287              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    2288            6 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: achint, acint, bchint, bcint
    2289              :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
    2290           24 :       TYPE(block_p_type), DIMENSION(3)                   :: blocks_rv
    2291            6 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
    2292              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    2293            6 :       TYPE(sap_int_type), DIMENSION(:), POINTER          :: sap_int
    2294              : 
    2295              : !$    INTEGER(kind=omp_lock_kind), &
    2296            6 : !$       ALLOCATABLE, DIMENSION(:) :: locks
    2297              : !$    INTEGER                                            :: lock_num, hash
    2298              : !$    INTEGER, PARAMETER                                 :: nlock = 501
    2299              : 
    2300            6 :       ppnl_present = ASSOCIATED(sap_ppnl)
    2301            6 :       IF (.NOT. ppnl_present) RETURN
    2302              : 
    2303            6 :       CALL timeset(routineN, handle)
    2304            6 :       nkind = SIZE(qs_kind_set)
    2305            6 :       natom = SIZE(particle_set)
    2306              : 
    2307              :       ! sap_int needs to be shared as multiple threads need to access this
    2308            6 :       NULLIFY (sap_int)
    2309           42 :       ALLOCATE (sap_int(nkind*nkind))
    2310           30 :       DO i = 1, nkind*nkind
    2311           24 :          NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
    2312           30 :          sap_int(i)%nalist = 0
    2313              :       END DO
    2314              : 
    2315              :       ! "nder" in moment_mode is "order"
    2316            6 :       CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=1, moment_mode=.TRUE., pseudoatom=pseudoatom)
    2317              : 
    2318              :       ! *** Set up a sorting index
    2319            6 :       CALL sap_sort(sap_int)
    2320              : 
    2321           30 :       ALLOCATE (basis_set(nkind))
    2322           18 :       DO ikind = 1, nkind
    2323           12 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
    2324           18 :          IF (ASSOCIATED(orb_basis_set)) THEN
    2325           12 :             basis_set(ikind)%gto_basis_set => orb_basis_set
    2326              :          ELSE
    2327            0 :             NULLIFY (basis_set(ikind)%gto_basis_set)
    2328              :          END IF
    2329              :       END DO
    2330              : 
    2331              :       ! *** All integrals needed have been calculated and stored in sap_int
    2332              :       ! *** We now calculate the commutator matrix elements
    2333              : 
    2334              : !$OMP PARALLEL &
    2335              : !$OMP DEFAULT (NONE) &
    2336              : !$OMP SHARED  (basis_set, matrix_rv, &
    2337              : !$OMP          sap_int, nkind, eps_ppnl, locks, sab_all) &
    2338              : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
    2339              : !$OMP          iab, irow, icol, blocks_rv, lock_num, &
    2340              : !$OMP          found, iac, ibc, alist_ac, alist_bc, &
    2341              : !$OMP          na, np, nb, kkind, kac, kbc, i, &
    2342            6 : !$OMP          hash, natom, acint, bcint, achint, bchint)
    2343              : 
    2344              : !$OMP SINGLE
    2345              : !$    ALLOCATE (locks(nlock))
    2346              : !$OMP END SINGLE
    2347              : 
    2348              : !$OMP DO
    2349              : !$    DO lock_num = 1, nlock
    2350              : !$       call omp_init_lock(locks(lock_num))
    2351              : !$    END DO
    2352              : !$OMP END DO
    2353              : 
    2354              : !$OMP DO SCHEDULE(GUIDED)
    2355              : 
    2356              :       DO slot = 1, sab_all(1)%nl_size
    2357              : 
    2358              :          ikind = sab_all(1)%nlist_task(slot)%ikind
    2359              :          jkind = sab_all(1)%nlist_task(slot)%jkind
    2360              :          iatom = sab_all(1)%nlist_task(slot)%iatom
    2361              :          jatom = sab_all(1)%nlist_task(slot)%jatom
    2362              :          cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
    2363              :          rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
    2364              : 
    2365              :          IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
    2366              :          IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
    2367              :          iab = ikind + nkind*(jkind - 1)
    2368              : 
    2369              :          ! *** Create matrix blocks for a new matrix block column ***
    2370              :          irow = iatom
    2371              :          icol = jatom
    2372              :          DO i = 1, 3
    2373              :             CALL dbcsr_get_block_p(matrix_rv(i)%matrix, irow, icol, blocks_rv(i)%block, found)
    2374              :             CPASSERT(found)
    2375              :          END DO
    2376              : 
    2377              :          ! loop over all kinds for projector atom
    2378              :          DO kkind = 1, nkind
    2379              :             iac = ikind + nkind*(kkind - 1)
    2380              :             ibc = jkind + nkind*(kkind - 1)
    2381              :             IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
    2382              :             IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
    2383              :             CALL get_alist(sap_int(iac), alist_ac, iatom)
    2384              :             CALL get_alist(sap_int(ibc), alist_bc, jatom)
    2385              : 
    2386              :             IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
    2387              :             IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
    2388              :             DO kac = 1, alist_ac%nclist
    2389              :                DO kbc = 1, alist_bc%nclist
    2390              :                   IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
    2391              :                   IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
    2392              :                      IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
    2393              :                      acint => alist_ac%clist(kac)%acint
    2394              :                      bcint => alist_bc%clist(kbc)%acint
    2395              :                      achint => alist_ac%clist(kac)%achint
    2396              :                      bchint => alist_bc%clist(kbc)%achint
    2397              :                      na = SIZE(acint, 1)
    2398              :                      np = SIZE(acint, 2)
    2399              :                      nb = SIZE(bcint, 1)
    2400              : !$                   hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
    2401              : !$                   CALL omp_set_lock(locks(hash))
    2402              :                      ! Vnl*r
    2403              :                      blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
    2404              :                                                       MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
    2405              :                      blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
    2406              :                                                       MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
    2407              :                      blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
    2408              :                                                       MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
    2409              : 
    2410              :                      ! r*Vnl
    2411              :                      blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) - &
    2412              :                                                       MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2413              :                      blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) - &
    2414              :                                                       MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2415              :                      blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) - &
    2416              :                                                       MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
    2417              : 
    2418              : !$                   CALL omp_unset_lock(locks(hash))
    2419              :                      EXIT ! We have found a match and there can be only one single match
    2420              :                   END IF
    2421              :                END DO
    2422              :             END DO
    2423              :          END DO
    2424              :          DO i = 1, 3
    2425              :             NULLIFY (blocks_rv(i)%block)
    2426              :          END DO
    2427              :       END DO
    2428              : 
    2429              : !$OMP DO
    2430              : !$    DO lock_num = 1, nlock
    2431              : !$       call omp_destroy_lock(locks(lock_num))
    2432              : !$    END DO
    2433              : !$OMP END DO
    2434              : 
    2435              : !$OMP SINGLE
    2436              : !$    DEALLOCATE (locks)
    2437              : !$OMP END SINGLE NOWAIT
    2438              : 
    2439              : !$OMP END PARALLEL
    2440              : 
    2441            6 :       CALL release_sap_int(sap_int)
    2442              : 
    2443            6 :       DEALLOCATE (basis_set)
    2444              : 
    2445            6 :       CALL timestop(handle)
    2446              : 
    2447           18 :    END SUBROUTINE build_drpnl_matrix
    2448              : 
    2449              : ! **************************************************************************************************
    2450              : !> \brief   Adapted from comab_opr. Calculate the product O*r or r*O from the integrals [a|O|b].
    2451              : !>          We assume that on input all integrals [a+1|O|b+1] are available.
    2452              : !> \param la_max ...
    2453              : !> \param npgfa ...
    2454              : !> \param rpgfa ...
    2455              : !> \param la_min ...
    2456              : !> \param lb_max ...
    2457              : !> \param npgfb ...
    2458              : !> \param rpgfb ...
    2459              : !> \param lb_min ...
    2460              : !> \param dab ...
    2461              : !> \param ab ...
    2462              : !> \param comabr ...
    2463              : !>
    2464              : !> \param ra ...
    2465              : !> \param rb ...
    2466              : !> \param direction_Or ...
    2467              : !> \par Literature
    2468              : !>          S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
    2469              : !> \par Parameters
    2470              : !>      - ax,ay,az  : Angular momentum index numbers of orbital a.
    2471              : !>      - bx,by,bz  : Angular momentum index numbers of orbital b.
    2472              : !>      - coset     : Cartesian orbital set pointer.
    2473              : !>      - l{a,b}    : Angular momentum quantum number of shell a or b.
    2474              : !>      - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
    2475              : !>      - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
    2476              : !>      - ncoset    : Number of orbitals in a Cartesian orbital set.
    2477              : !>      - npgf{a,b} : Degree of contraction of shell a or b.
    2478              : !>      - rab       : Distance vector between the atomic centers a and b.
    2479              : !>      - rab2      : Square of the distance between the atomic centers a and b.
    2480              : !>      - rac       : Distance vector between the atomic centers a and c.
    2481              : !>      - rac2      : Square of the distance between the atomic centers a and c.
    2482              : !>      - rbc       : Distance vector between the atomic centers b and c.
    2483              : !>      - rbc2      : Square of the distance between the atomic centers b and c.
    2484              : !>      - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
    2485              : !>      - zet{a,b}  : Exponents of the Gaussian-type functions a or b.
    2486              : !>      - zetp      : Reciprocal of the sum of the exponents of orbital a and b.
    2487              : !>
    2488              : !> \author  Tomas Zimmermann
    2489              : ! **************************************************************************************************
    2490         2052 :    SUBROUTINE ab_opr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
    2491         4104 :                      dab, ab, comabr, ra, rb, direction_Or)
    2492              :       INTEGER, INTENT(IN)                                :: la_max, npgfa
    2493              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa
    2494              :       INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
    2495              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb
    2496              :       INTEGER, INTENT(IN)                                :: lb_min
    2497              :       REAL(KIND=dp), INTENT(IN)                          :: dab
    2498              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: ab
    2499              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: comabr
    2500              :       REAL(KIND=dp), DIMENSION(1:3), INTENT(IN)          :: ra, rb
    2501              :       LOGICAL                                            :: direction_Or
    2502              : 
    2503              :       INTEGER                                            :: ax, ay, az, bx, by, bz, coa, coap, &
    2504              :                                                             coapx, coapy, coapz, cob, cobp, cobpx, &
    2505              :                                                             cobpy, cobpz, ipgf, jpgf, la, lb, na, &
    2506              :                                                             nap, nb, nbp, ofa, ofb
    2507              : 
    2508     82941840 :       comabr = 0.0_dp
    2509              : 
    2510         2052 :       ofa = ncoset(la_min - 1)
    2511         2052 :       ofb = ncoset(lb_min - 1)
    2512              : 
    2513         2052 :       na = 0
    2514         2052 :       nap = 0
    2515        10260 :       DO ipgf = 1, npgfa
    2516         8208 :          nb = 0
    2517         8208 :          nbp = 0
    2518        41040 :          DO jpgf = 1, npgfb
    2519        32832 :             IF (rpgfa(ipgf) + rpgfb(jpgf) > dab) THEN
    2520        77596 :                DO la = la_min, la_max
    2521       135546 :                   DO ax = 0, la
    2522       173850 :                      DO ay = 0, la - ax
    2523        70718 :                         az = la - ax - ay
    2524        70718 :                         coa = na + coset(ax, ay, az) - ofa
    2525        70718 :                         coap = nap + coset(ax, ay, az) - ofa
    2526        70718 :                         coapx = nap + coset(ax + 1, ay, az) - ofa
    2527        70718 :                         coapy = nap + coset(ax, ay + 1, az) - ofa
    2528        70718 :                         coapz = nap + coset(ax, ay, az + 1) - ofa
    2529       221274 :                         DO lb = lb_min, lb_max
    2530       277818 :                            DO bx = 0, lb
    2531       343482 :                               DO by = 0, lb - bx
    2532       136382 :                                  bz = lb - bx - by
    2533       136382 :                                  cob = nb + coset(bx, by, bz) - ofb
    2534       136382 :                                  cobp = nbp + coset(bx, by, bz) - ofb
    2535       136382 :                                  cobpx = nbp + coset(bx + 1, by, bz) - ofb
    2536       136382 :                                  cobpy = nbp + coset(bx, by + 1, bz) - ofb
    2537       136382 :                                  cobpz = nbp + coset(bx, by, bz + 1) - ofb
    2538       250876 :                                  IF (direction_Or) THEN
    2539              :                                     ! [a|O * x|b] = [a|O|b(x+1)] + [a|O|b] * X_b
    2540              :                                     !             = [a|O * (x - X_b)|b] + [a|O|b] * X_b
    2541              :                                     ! So the second term makes sure that we actually calculate
    2542              :                                     !  <O*r> and not <O*(r-R)>
    2543        19912 :                                     comabr(coa, cob, 1) = ab(coap, cobpx) + ab(coap, cobp)*rb(1)
    2544        19912 :                                     comabr(coa, cob, 2) = ab(coap, cobpy) + ab(coap, cobp)*rb(2)
    2545        19912 :                                     comabr(coa, cob, 3) = ab(coap, cobpz) + ab(coap, cobp)*rb(3)
    2546              :                                  ELSE
    2547       116470 :                                     comabr(coa, cob, 1) = ab(coapx, cobp) + ab(coap, cobp)*ra(1)
    2548       116470 :                                     comabr(coa, cob, 2) = ab(coapy, cobp) + ab(coap, cobp)*ra(2)
    2549       116470 :                                     comabr(coa, cob, 3) = ab(coapz, cobp) + ab(coap, cobp)*ra(3)
    2550              :                                  END IF
    2551              :                               END DO
    2552              :                            END DO
    2553              :                         END DO
    2554              :                      END DO
    2555              :                   END DO
    2556              :                END DO
    2557              :             END IF
    2558        32832 :             nb = nb + ncoset(lb_max) - ofb
    2559        41040 :             nbp = nbp + ncoset(lb_max + 1) - ofb
    2560              :          END DO
    2561         8208 :          na = na + ncoset(la_max) - ofa
    2562        10260 :          nap = nap + ncoset(la_max + 1) - ofa
    2563              :       END DO
    2564              : 
    2565         2052 :    END SUBROUTINE ab_opr
    2566              : 
    2567              : ! **************************************************************************************************
    2568              : !> \brief Apply the operator \delta_\mu^\lambda to zero out all elements of the matrix
    2569              : !>         which don't fulfill the condition.
    2570              : !> \param matrix ...
    2571              : !> \param qs_kind_set ...
    2572              : !> \param basis_type ...
    2573              : !> \param sab_nl ...
    2574              : !> \param lambda ...
    2575              : !> \param direction_Or ...
    2576              : !> \author Edward Ditler
    2577              : ! **************************************************************************************************
    2578            0 :    SUBROUTINE hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
    2579              : 
    2580              :       TYPE(dbcsr_type)                                   :: matrix
    2581              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2582              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
    2583              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2584              :          POINTER                                         :: sab_nl
    2585              :       INTEGER, INTENT(IN)                                :: lambda
    2586              :       LOGICAL                                            :: direction_Or
    2587              : 
    2588              :       CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_1d'
    2589              : 
    2590              :       INTEGER                                            :: handle, iatom, icol, ikind, irow, jatom, &
    2591              :                                                             jkind, ldsab, mepos, nkind, nseta, &
    2592              :                                                             nsetb, nthread
    2593              :       INTEGER, DIMENSION(3)                              :: cell
    2594            0 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    2595            0 :                                                             npgfb, nsgfa, nsgfb
    2596            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    2597              :       LOGICAL                                            :: do_symmetric, found
    2598              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    2599            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
    2600            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: k_block, rpgfa, rpgfb, scon_a, scon_b, &
    2601            0 :                                                             sphi_a, sphi_b, zeta, zetb
    2602            0 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
    2603              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
    2604              :       TYPE(neighbor_list_iterator_p_type), &
    2605            0 :          DIMENSION(:), POINTER                           :: nl_iterator
    2606              : 
    2607            0 :       CALL timeset(routineN, handle)
    2608              : 
    2609            0 :       nkind = SIZE(qs_kind_set)
    2610              : 
    2611              :       ! check for symmetry
    2612            0 :       CPASSERT(SIZE(sab_nl) > 0)
    2613            0 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
    2614              : 
    2615              :       ! prepare basis set
    2616            0 :       ALLOCATE (basis_set_list(nkind))
    2617            0 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
    2618              : 
    2619              :       ! *** Allocate work storage ***
    2620            0 :       ldsab = get_memory_usage(qs_kind_set, basis_type)
    2621              : 
    2622              :       nthread = 1
    2623            0 : !$    nthread = omp_get_max_threads()
    2624              :       ! Iterate of neighbor list
    2625            0 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
    2626              : 
    2627              : !$OMP PARALLEL DEFAULT(NONE) &
    2628              : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
    2629              : !$OMP SHARED (ncoset,matrix,basis_set_list) &
    2630              : !$OMP SHARED (direction_or, lambda) &
    2631              : !$OMP PRIVATE (k_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
    2632              : !$OMP PRIVATE (basis_set_a,basis_set_b) &
    2633              : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
    2634              : !$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
    2635            0 : !$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found, sphi_a, sphi_b)
    2636              : 
    2637              :       mepos = 0
    2638              : !$    mepos = omp_get_thread_num()
    2639              : 
    2640              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
    2641              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
    2642              :                                 iatom=iatom, jatom=jatom, r=rab, cell=cell)
    2643              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
    2644              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    2645              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
    2646              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    2647              :          ! basis ikind
    2648              :          first_sgfa => basis_set_a%first_sgf
    2649              :          la_max => basis_set_a%lmax
    2650              :          la_min => basis_set_a%lmin
    2651              :          npgfa => basis_set_a%npgf
    2652              :          nsgfa => basis_set_a%nsgf_set
    2653              :          rpgfa => basis_set_a%pgf_radius
    2654              :          set_radius_a => basis_set_a%set_radius
    2655              :          sphi_a => basis_set_a%sphi
    2656              :          zeta => basis_set_a%zet
    2657              :          scon_a => basis_set_a%scon
    2658              :          ! basis jkind
    2659              :          first_sgfb => basis_set_b%first_sgf
    2660              :          lb_max => basis_set_b%lmax
    2661              :          lb_min => basis_set_b%lmin
    2662              :          npgfb => basis_set_b%npgf
    2663              :          nsgfb => basis_set_b%nsgf_set
    2664              :          rpgfb => basis_set_b%pgf_radius
    2665              :          set_radius_b => basis_set_b%set_radius
    2666              :          sphi_b => basis_set_b%sphi
    2667              :          zetb => basis_set_b%zet
    2668              :          scon_b => basis_set_b%scon
    2669              : 
    2670              :          nseta = basis_set_a%nset
    2671              :          nsetb = basis_set_b%nset
    2672              : 
    2673              :          IF (do_symmetric) THEN
    2674              :             IF (iatom <= jatom) THEN
    2675              :                irow = iatom
    2676              :                icol = jatom
    2677              :             ELSE
    2678              :                irow = jatom
    2679              :                icol = iatom
    2680              :             END IF
    2681              :          ELSE
    2682              :             irow = iatom
    2683              :             icol = jatom
    2684              :          END IF
    2685              : 
    2686              :          NULLIFY (k_block)
    2687              :          CALL dbcsr_get_block_p(matrix, irow, icol, k_block, found)
    2688              :          CPASSERT(found)
    2689              : 
    2690              :          IF (direction_Or) THEN
    2691              :             IF (jatom /= lambda) k_block(:, :) = 0._dp
    2692              :          ELSE IF (.NOT. direction_Or) THEN
    2693              :             IF (iatom /= lambda) k_block(:, :) = 0._dp
    2694              :          END IF
    2695              :       END DO
    2696              : !$OMP END PARALLEL
    2697            0 :       CALL neighbor_list_iterator_release(nl_iterator)
    2698              : 
    2699              :       ! Release work storage
    2700            0 :       DEALLOCATE (basis_set_list)
    2701              : 
    2702            0 :       CALL timestop(handle)
    2703              : 
    2704            0 :    END SUBROUTINE hr_mult_by_delta_1d
    2705              : 
    2706              : ! **************************************************************************************************
    2707              : !> \brief Apply the operator \delta_\mu^\lambda to zero out all elements of the matrix
    2708              : !>         which don't fulfill the condition.
    2709              : !>        Operates on matrix_hr(1:3) instead of a single matrix
    2710              : !> \param matrix_hr ...
    2711              : !> \param qs_kind_set ...
    2712              : !> \param basis_type ...
    2713              : !> \param sab_nl ...
    2714              : !> \param deltaR ...
    2715              : !> \param direction_Or ...
    2716              : !> \author Edward Ditler
    2717              : ! **************************************************************************************************
    2718           18 :    SUBROUTINE hr_mult_by_delta_3d(matrix_hr, qs_kind_set, basis_type, sab_nl, deltaR, direction_Or)
    2719              : 
    2720              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hr
    2721              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2722              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
    2723              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2724              :          POINTER                                         :: sab_nl
    2725              :       REAL(KIND=dp), DIMENSION(:, :)                     :: deltaR
    2726              :       LOGICAL                                            :: direction_Or
    2727              : 
    2728              :       CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_3d'
    2729              : 
    2730              :       INTEGER                                            :: handle, iatom, icol, ikind, ir, irow, &
    2731              :                                                             jatom, jkind, ldsab, mepos, nkind, &
    2732              :                                                             nseta, nsetb, nthread
    2733              :       INTEGER, DIMENSION(3)                              :: cell
    2734           18 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    2735           18 :                                                             npgfb, nsgfa, nsgfb
    2736           18 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    2737              :       LOGICAL                                            :: do_symmetric, found
    2738              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    2739           18 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
    2740           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: kx_block, ky_block, kz_block, rpgfa, &
    2741           18 :                                                             rpgfb, scon_a, scon_b, sphi_a, sphi_b, &
    2742           18 :                                                             zeta, zetb
    2743           18 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
    2744              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
    2745              :       TYPE(neighbor_list_iterator_p_type), &
    2746           18 :          DIMENSION(:), POINTER                           :: nl_iterator
    2747              : 
    2748           18 :       CALL timeset(routineN, handle)
    2749              : 
    2750           18 :       nkind = SIZE(qs_kind_set)
    2751              : 
    2752              :       ! check for symmetry
    2753           18 :       CPASSERT(SIZE(sab_nl) > 0)
    2754           18 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
    2755              : 
    2756              :       ! prepare basis set
    2757           90 :       ALLOCATE (basis_set_list(nkind))
    2758           18 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
    2759              : 
    2760              :       ! *** Allocate work storage ***
    2761           18 :       ldsab = get_memory_usage(qs_kind_set, basis_type)
    2762              : 
    2763              :       nthread = 1
    2764           18 : !$    nthread = omp_get_max_threads()
    2765              :       ! Iterate of neighbor list
    2766           18 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
    2767              : 
    2768              : !$OMP PARALLEL DEFAULT(NONE) &
    2769              : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
    2770              : !$OMP SHARED (ncoset,matrix_hr,basis_set_list) &
    2771              : !$OMP SHARED (direction_or, deltar) &
    2772              : !$OMP PRIVATE (kx_block,ky_block,kz_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
    2773              : !$OMP PRIVATE (basis_set_a,basis_set_b) &
    2774              : !$OMP PRIVATE (nseta, nsetb) &
    2775              : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, scon_a) &
    2776              : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, sphi_b, zetb, scon_b) &
    2777           18 : !$OMP PRIVATE (irow, icol, found)
    2778              : 
    2779              :       mepos = 0
    2780              : !$    mepos = omp_get_thread_num()
    2781              : 
    2782              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
    2783              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
    2784              :                                 iatom=iatom, jatom=jatom, r=rab, cell=cell)
    2785              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
    2786              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    2787              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
    2788              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    2789              :          ! basis ikind
    2790              :          first_sgfa => basis_set_a%first_sgf
    2791              :          la_max => basis_set_a%lmax
    2792              :          la_min => basis_set_a%lmin
    2793              :          npgfa => basis_set_a%npgf
    2794              :          nsgfa => basis_set_a%nsgf_set
    2795              :          rpgfa => basis_set_a%pgf_radius
    2796              :          set_radius_a => basis_set_a%set_radius
    2797              :          sphi_a => basis_set_a%sphi
    2798              :          zeta => basis_set_a%zet
    2799              :          scon_a => basis_set_a%scon
    2800              :          ! basis jkind
    2801              :          first_sgfb => basis_set_b%first_sgf
    2802              :          lb_max => basis_set_b%lmax
    2803              :          lb_min => basis_set_b%lmin
    2804              :          npgfb => basis_set_b%npgf
    2805              :          nsgfb => basis_set_b%nsgf_set
    2806              :          rpgfb => basis_set_b%pgf_radius
    2807              :          set_radius_b => basis_set_b%set_radius
    2808              :          sphi_b => basis_set_b%sphi
    2809              :          zetb => basis_set_b%zet
    2810              :          scon_b => basis_set_b%scon
    2811              : 
    2812              :          nseta = basis_set_a%nset
    2813              :          nsetb = basis_set_b%nset
    2814              : 
    2815              :          IF (do_symmetric) THEN
    2816              :             IF (iatom <= jatom) THEN
    2817              :                irow = iatom
    2818              :                icol = jatom
    2819              :             ELSE
    2820              :                irow = jatom
    2821              :                icol = iatom
    2822              :             END IF
    2823              :          ELSE
    2824              :             irow = iatom
    2825              :             icol = jatom
    2826              :          END IF
    2827              : 
    2828              :          NULLIFY (kx_block, ky_block, kz_block)
    2829              :          CALL dbcsr_get_block_p(matrix_hr(1)%matrix, irow, icol, kx_block, found)
    2830              :          CPASSERT(found)
    2831              :          CALL dbcsr_get_block_p(matrix_hr(2)%matrix, irow, icol, ky_block, found)
    2832              :          CPASSERT(found)
    2833              :          CALL dbcsr_get_block_p(matrix_hr(3)%matrix, irow, icol, kz_block, found)
    2834              :          CPASSERT(found)
    2835              : 
    2836              :          IF (direction_Or) THEN
    2837              :             DO ir = 1, 3
    2838              : !$OMP CRITICAL(blockadd)
    2839              :                SELECT CASE (ir)
    2840              :                CASE (1)
    2841              :                   kx_block(:, :) = kx_block(:, :)*deltaR(ir, jatom)
    2842              :                CASE (2)
    2843              :                   ky_block(:, :) = ky_block(:, :)*deltaR(ir, jatom)
    2844              :                CASE (3)
    2845              :                   kz_block(:, :) = kz_block(:, :)*deltaR(ir, jatom)
    2846              :                END SELECT
    2847              : !$OMP END CRITICAL(blockadd)
    2848              :             END DO
    2849              :          ELSE
    2850              :             DO ir = 1, 3
    2851              : !$OMP CRITICAL(blockadd)
    2852              :                SELECT CASE (ir)
    2853              :                CASE (1)
    2854              :                   kx_block(:, :) = kx_block(:, :)*deltaR(ir, iatom)
    2855              :                CASE (2)
    2856              :                   ky_block(:, :) = ky_block(:, :)*deltaR(ir, iatom)
    2857              :                CASE (3)
    2858              :                   kz_block(:, :) = kz_block(:, :)*deltaR(ir, iatom)
    2859              :                END SELECT
    2860              : !$OMP END CRITICAL(blockadd)
    2861              :             END DO
    2862              :          END IF
    2863              :       END DO
    2864              : !$OMP END PARALLEL
    2865           18 :       CALL neighbor_list_iterator_release(nl_iterator)
    2866              : 
    2867              :       ! Release work storage
    2868           18 :       DEALLOCATE (basis_set_list)
    2869              : 
    2870           18 :       CALL timestop(handle)
    2871              : 
    2872           36 :    END SUBROUTINE hr_mult_by_delta_3d
    2873              : 
    2874         4104 : END MODULE qs_vcd_ao
    2875              : 
        

Generated by: LCOV version 2.0-1