LCOV - code coverage report
Current view: top level - src - qs_o3c_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 53 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 3 0

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : ! **************************************************************************************************
       8              : !> \brief Methods used with 3-center overlap type integrals containers
       9              : !> \par History
      10              : !>      - none
      11              : !>      - 11.2018 fixed OMP race condition in contract3_o3c routine (A.Bussy)
      12              : ! **************************************************************************************************
      13              : MODULE qs_o3c_methods
      14              :    USE ai_contraction_sphi,             ONLY: abc_contract
      15              :    USE ai_overlap3,                     ONLY: overlap3
      16              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      17              :                                               gto_basis_set_type
      18              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
      19              :                                               dbcsr_p_type,&
      20              :                                               dbcsr_type
      21              :    USE kinds,                           ONLY: dp
      22              :    USE orbital_pointers,                ONLY: ncoset
      23              :    USE qs_o3c_types,                    ONLY: &
      24              :         get_o3c_container, get_o3c_iterator_info, get_o3c_vec, o3c_container_type, o3c_iterate, &
      25              :         o3c_iterator_create, o3c_iterator_release, o3c_iterator_type, o3c_vec_type, &
      26              :         set_o3c_container
      27              : 
      28              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
      29              : #include "./base/base_uses.f90"
      30              : 
      31              :    IMPLICIT NONE
      32              : 
      33              :    PRIVATE
      34              : 
      35              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_o3c_methods'
      36              : 
      37              :    PUBLIC :: calculate_o3c_integrals, contract12_o3c, contract3_o3c
      38              : 
      39              : CONTAINS
      40              : 
      41              : ! **************************************************************************************************
      42              : !> \brief ...
      43              : !> \param o3c ...
      44              : !> \param calculate_forces ...
      45              : !> \param matrix_p ...
      46              : ! **************************************************************************************************
      47            0 :    SUBROUTINE calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
      48              :       TYPE(o3c_container_type), POINTER                  :: o3c
      49              :       LOGICAL, INTENT(IN), OPTIONAL                      :: calculate_forces
      50              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
      51              :          POINTER                                         :: matrix_p
      52              : 
      53              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_o3c_integrals'
      54              : 
      55              :       INTEGER :: egfa, egfb, egfc, handle, i, iatom, icol, ikind, irow, iset, ispin, j, jatom, &
      56              :          jkind, jset, katom, kkind, kset, mepos, ncoa, ncob, ncoc, ni, nj, nk, nseta, nsetb, &
      57              :          nsetc, nspin, nthread, sgfa, sgfb, sgfc
      58            0 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, lc_max, &
      59            0 :                                                             lc_min, npgfa, npgfb, npgfc, nsgfa, &
      60            0 :                                                             nsgfb, nsgfc
      61            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfc
      62              :       LOGICAL                                            :: do_force, found, trans
      63              :       REAL(KIND=dp)                                      :: dij, dik, djk, fpre
      64            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: pmat
      65            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: sabc, sabc_contr
      66            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: iabdc, iadbc, idabc, sabdc, sdabc
      67              :       REAL(KIND=dp), DIMENSION(3)                        :: rij, rik, rjk
      68            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_c
      69            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fi, fj, fk, pblock, rpgfa, rpgfb, rpgfc, &
      70            0 :                                                             sphi_a, sphi_b, sphi_c, tvec, zeta, &
      71            0 :                                                             zetb, zetc
      72            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: iabc
      73            0 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list_a, basis_set_list_b, &
      74            0 :                                                             basis_set_list_c
      75              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b, basis_set_c
      76              :       TYPE(o3c_iterator_type)                            :: o3c_iterator
      77              : 
      78            0 :       CALL timeset(routineN, handle)
      79              : 
      80            0 :       do_force = .FALSE.
      81            0 :       IF (PRESENT(calculate_forces)) do_force = calculate_forces
      82            0 :       CALL get_o3c_container(o3c, nspin=nspin)
      83              : 
      84              :       ! basis sets
      85              :       CALL get_o3c_container(o3c, basis_set_list_a=basis_set_list_a, &
      86            0 :                              basis_set_list_b=basis_set_list_b, basis_set_list_c=basis_set_list_c)
      87              : 
      88              :       nthread = 1
      89            0 : !$    nthread = omp_get_max_threads()
      90            0 :       CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
      91              : 
      92              : !$OMP PARALLEL DEFAULT(NONE) &
      93              : !$OMP SHARED (nthread,o3c_iterator,ncoset,nspin,basis_set_list_a,basis_set_list_b,&
      94              : !$OMP         basis_set_list_c,do_force,matrix_p)&
      95              : !$OMP PRIVATE (mepos,ikind,jkind,kkind,basis_set_a,basis_set_b,basis_set_c,rij,rik,rjk,&
      96              : !$OMP          first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,rpgfa,set_radius_a,sphi_a,zeta,&
      97              : !$OMP          first_sgfb,lb_max,lb_min,npgfb,nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,&
      98              : !$OMP          first_sgfc,lc_max,lc_min,npgfc,nsetc,nsgfc,rpgfc,set_radius_c,sphi_c,zetc,&
      99              : !$OMP          iset,jset,kset,dij,dik,djk,ni,nj,nk,iabc,idabc,iadbc,iabdc,tvec,fi,fj,fk,ncoa,&
     100              : !$OMP          ncob,ncoc,sabc,sabc_contr,sdabc,sabdc,sgfa,sgfb,sgfc,egfa,egfb,egfc,i,j,&
     101            0 : !$OMP          pblock,pmat,ispin,iatom,jatom,katom,irow,icol,found,trans,fpre)
     102              : 
     103              :       mepos = 0
     104              : !$    mepos = omp_get_thread_num()
     105              : 
     106              :       DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
     107              :          CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, &
     108              :                                     ikind=ikind, jkind=jkind, kkind=kkind, rij=rij, rik=rik, &
     109              :                                     integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
     110              :          CPASSERT(.NOT. ASSOCIATED(iabc))
     111              :          CPASSERT(.NOT. ASSOCIATED(tvec))
     112              :          CPASSERT(.NOT. ASSOCIATED(fi))
     113              :          CPASSERT(.NOT. ASSOCIATED(fj))
     114              :          CPASSERT(.NOT. ASSOCIATED(fk))
     115              :          ! basis
     116              :          basis_set_a => basis_set_list_a(ikind)%gto_basis_set
     117              :          basis_set_b => basis_set_list_b(jkind)%gto_basis_set
     118              :          basis_set_c => basis_set_list_c(kkind)%gto_basis_set
     119              :          ! center A
     120              :          first_sgfa => basis_set_a%first_sgf
     121              :          la_max => basis_set_a%lmax
     122              :          la_min => basis_set_a%lmin
     123              :          npgfa => basis_set_a%npgf
     124              :          nseta = basis_set_a%nset
     125              :          nsgfa => basis_set_a%nsgf_set
     126              :          rpgfa => basis_set_a%pgf_radius
     127              :          set_radius_a => basis_set_a%set_radius
     128              :          sphi_a => basis_set_a%sphi
     129              :          zeta => basis_set_a%zet
     130              :          ! center B
     131              :          first_sgfb => basis_set_b%first_sgf
     132              :          lb_max => basis_set_b%lmax
     133              :          lb_min => basis_set_b%lmin
     134              :          npgfb => basis_set_b%npgf
     135              :          nsetb = basis_set_b%nset
     136              :          nsgfb => basis_set_b%nsgf_set
     137              :          rpgfb => basis_set_b%pgf_radius
     138              :          set_radius_b => basis_set_b%set_radius
     139              :          sphi_b => basis_set_b%sphi
     140              :          zetb => basis_set_b%zet
     141              :          ! center C (RI)
     142              :          first_sgfc => basis_set_c%first_sgf
     143              :          lc_max => basis_set_c%lmax
     144              :          lc_min => basis_set_c%lmin
     145              :          npgfc => basis_set_c%npgf
     146              :          nsetc = basis_set_c%nset
     147              :          nsgfc => basis_set_c%nsgf_set
     148              :          rpgfc => basis_set_c%pgf_radius
     149              :          set_radius_c => basis_set_c%set_radius
     150              :          sphi_c => basis_set_c%sphi
     151              :          zetc => basis_set_c%zet
     152              : 
     153              :          ni = SUM(nsgfa)
     154              :          nj = SUM(nsgfb)
     155              :          nk = SUM(nsgfc)
     156              : 
     157              :          ALLOCATE (iabc(ni, nj, nk))
     158              :          iabc(:, :, :) = 0.0_dp
     159              :          IF (do_force) THEN
     160              :             ALLOCATE (fi(nk, 3), fj(nk, 3), fk(nk, 3))
     161              :             fi(:, :) = 0.0_dp
     162              :             fj(:, :) = 0.0_dp
     163              :             fk(:, :) = 0.0_dp
     164              :             ALLOCATE (idabc(ni, nj, nk, 3))
     165              :             idabc(:, :, :, :) = 0.0_dp
     166              :             ALLOCATE (iadbc(ni, nj, nk, 3))
     167              :             iadbc(:, :, :, :) = 0.0_dp
     168              :             ALLOCATE (iabdc(ni, nj, nk, 3))
     169              :             iabdc(:, :, :, :) = 0.0_dp
     170              :          ELSE
     171              :             NULLIFY (fi, fj, fk)
     172              :          END IF
     173              :          ALLOCATE (tvec(nk, nspin))
     174              :          tvec(:, :) = 0.0_dp
     175              : 
     176              :          rjk(1:3) = rik(1:3) - rij(1:3)
     177              :          dij = NORM2(rij)
     178              :          dik = NORM2(rik)
     179              :          djk = NORM2(rjk)
     180              : 
     181              :          DO iset = 1, nseta
     182              :             DO jset = 1, nsetb
     183              :                IF (set_radius_a(iset) + set_radius_b(jset) < dij) CYCLE
     184              :                DO kset = 1, nsetc
     185              :                   IF (set_radius_a(iset) + set_radius_c(kset) < dik) CYCLE
     186              :                   IF (set_radius_b(jset) + set_radius_c(kset) < djk) CYCLE
     187              : 
     188              :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     189              :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
     190              :                   ncoc = npgfc(kset)*ncoset(lc_max(kset))
     191              : 
     192              :                   sgfa = first_sgfa(1, iset)
     193              :                   sgfb = first_sgfb(1, jset)
     194              :                   sgfc = first_sgfc(1, kset)
     195              : 
     196              :                   egfa = sgfa + nsgfa(iset) - 1
     197              :                   egfb = sgfb + nsgfb(jset) - 1
     198              :                   egfc = sgfc + nsgfc(kset) - 1
     199              : 
     200              :                   IF (ncoa*ncob*ncoc > 0) THEN
     201              :                      ALLOCATE (sabc(ncoa, ncob, ncoc))
     202              :                      sabc(:, :, :) = 0.0_dp
     203              :                      IF (do_force) THEN
     204              :                         ALLOCATE (sdabc(ncoa, ncob, ncoc, 3))
     205              :                         sdabc(:, :, :, :) = 0.0_dp
     206              :                         ALLOCATE (sabdc(ncoa, ncob, ncoc, 3))
     207              :                         sabdc(:, :, :, :) = 0.0_dp
     208              :                         CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     209              :                                       lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     210              :                                       lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
     211              :                                       rij, dij, rik, dik, rjk, djk, sabc, sdabc, sabdc)
     212              :                      ELSE
     213              :                         CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     214              :                                       lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     215              :                                       lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
     216              :                                       rij, dij, rik, dik, rjk, djk, sabc)
     217              :                      END IF
     218              :                      ALLOCATE (sabc_contr(nsgfa(iset), nsgfb(jset), nsgfc(kset)))
     219              : 
     220              :                      CALL abc_contract(sabc_contr, sabc, &
     221              :                                        sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
     222              :                                        ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
     223              :                      iabc(sgfa:egfa, sgfb:egfb, sgfc:egfc) = &
     224              :                         sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
     225              :                      IF (do_force) THEN
     226              :                         DO i = 1, 3
     227              :                            CALL abc_contract(sabc_contr, sdabc(:, :, :, i), &
     228              :                                              sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
     229              :                                              ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
     230              :                            idabc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
     231              :                               sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
     232              :                            CALL abc_contract(sabc_contr, sabdc(:, :, :, i), &
     233              :                                              sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
     234              :                                              ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
     235              :                            iabdc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
     236              :                               sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
     237              :                         END DO
     238              :                      END IF
     239              : 
     240              :                      DEALLOCATE (sabc_contr)
     241              :                      DEALLOCATE (sabc)
     242              :                   END IF
     243              :                   IF (do_force) THEN
     244              :                      DEALLOCATE (sdabc, sabdc)
     245              :                   END IF
     246              :                END DO
     247              :             END DO
     248              :          END DO
     249              :          IF (do_force) THEN
     250              :             ! translational invariance
     251              :             iadbc(:, :, :, :) = -idabc(:, :, :, :) - iabdc(:, :, :, :)
     252              :             !
     253              :             ! get the atom indices
     254              :             CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, &
     255              :                                        iatom=iatom, jatom=jatom, katom=katom)
     256              :             !
     257              :             ! contract over i and j to get forces
     258              :             IF (iatom <= jatom) THEN
     259              :                irow = iatom
     260              :                icol = jatom
     261              :                trans = .FALSE.
     262              :             ELSE
     263              :                irow = jatom
     264              :                icol = iatom
     265              :                trans = .TRUE.
     266              :             END IF
     267              :             IF (iatom == jatom) THEN
     268              :                fpre = 1.0_dp
     269              :             ELSE
     270              :                fpre = 2.0_dp
     271              :             END IF
     272              :             ALLOCATE (pmat(ni, nj))
     273              :             pmat(:, :) = 0.0_dp
     274              :             DO ispin = 1, nspin
     275              :                CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
     276              :                                       row=irow, col=icol, BLOCK=pblock, found=found)
     277              :                IF (found) THEN
     278              :                   IF (trans) THEN
     279              :                      pmat(:, :) = pmat(:, :) + TRANSPOSE(pblock(:, :))
     280              :                   ELSE
     281              :                      pmat(:, :) = pmat(:, :) + pblock(:, :)
     282              :                   END IF
     283              :                END IF
     284              :             END DO
     285              :             DO i = 1, 3
     286              :                DO j = 1, nk
     287              :                   fi(j, i) = fpre*SUM(pmat(:, :)*idabc(:, :, j, i))
     288              :                   fj(j, i) = fpre*SUM(pmat(:, :)*iadbc(:, :, j, i))
     289              :                   fk(j, i) = fpre*SUM(pmat(:, :)*iabdc(:, :, j, i))
     290              :                END DO
     291              :             END DO
     292              :             DEALLOCATE (pmat)
     293              :             !
     294              :             DEALLOCATE (idabc, iadbc, iabdc)
     295              :          END IF
     296              :          !
     297              :          CALL set_o3c_container(o3c_iterator, mepos=mepos, &
     298              :                                 integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
     299              : 
     300              :       END DO
     301              : !$OMP END PARALLEL
     302            0 :       CALL o3c_iterator_release(o3c_iterator)
     303              : 
     304            0 :       CALL timestop(handle)
     305              : 
     306            0 :    END SUBROUTINE calculate_o3c_integrals
     307              : 
     308              : ! **************************************************************************************************
     309              : !> \brief Contraction of 3-tensor over indices 1 and 2 (assuming symmetry)
     310              : !>        t(k) = sum_ij (ijk)*p(ij)
     311              : !> \param o3c ...
     312              : !> \param matrix_p ...
     313              : ! **************************************************************************************************
     314            0 :    SUBROUTINE contract12_o3c(o3c, matrix_p)
     315              :       TYPE(o3c_container_type), POINTER                  :: o3c
     316              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p
     317              : 
     318              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract12_o3c'
     319              : 
     320              :       INTEGER                                            :: handle, iatom, icol, ik, irow, ispin, &
     321              :                                                             jatom, mepos, nk, nspin, nthread
     322              :       LOGICAL                                            :: found, ijsymmetric, trans
     323              :       REAL(KIND=dp)                                      :: fpre
     324            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pblock, tvec
     325            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: iabc
     326              :       TYPE(o3c_iterator_type)                            :: o3c_iterator
     327              : 
     328            0 :       CALL timeset(routineN, handle)
     329              : 
     330            0 :       nspin = SIZE(matrix_p, 1)
     331            0 :       CALL get_o3c_container(o3c, ijsymmetric=ijsymmetric)
     332            0 :       CPASSERT(ijsymmetric)
     333              : 
     334              :       nthread = 1
     335            0 : !$    nthread = omp_get_max_threads()
     336            0 :       CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
     337              : 
     338              : !$OMP PARALLEL DEFAULT(NONE) &
     339              : !$OMP SHARED (nthread,o3c_iterator,matrix_p,nspin)&
     340            0 : !$OMP PRIVATE (mepos,ispin,iatom,jatom,ik,nk,irow,icol,iabc,tvec,found,pblock,trans,fpre)
     341              : 
     342              :       mepos = 0
     343              : !$    mepos = omp_get_thread_num()
     344              : 
     345              :       DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
     346              :          CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, iatom=iatom, jatom=jatom, &
     347              :                                     integral=iabc, tvec=tvec)
     348              :          nk = SIZE(tvec, 1)
     349              : 
     350              :          IF (iatom <= jatom) THEN
     351              :             irow = iatom
     352              :             icol = jatom
     353              :             trans = .FALSE.
     354              :          ELSE
     355              :             irow = jatom
     356              :             icol = iatom
     357              :             trans = .TRUE.
     358              :          END IF
     359              :          IF (iatom == jatom) THEN
     360              :             fpre = 1.0_dp
     361              :          ELSE
     362              :             fpre = 2.0_dp
     363              :          END IF
     364              : 
     365              :          DO ispin = 1, nspin
     366              :             CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
     367              :                                    row=irow, col=icol, BLOCK=pblock, found=found)
     368              :             IF (found) THEN
     369              :                IF (trans) THEN
     370              :                   DO ik = 1, nk
     371              :                      tvec(ik, ispin) = fpre*SUM(TRANSPOSE(pblock(:, :))*iabc(:, :, ik))
     372              :                   END DO
     373              :                ELSE
     374              :                   DO ik = 1, nk
     375              :                      tvec(ik, ispin) = fpre*SUM(pblock(:, :)*iabc(:, :, ik))
     376              :                   END DO
     377              :                END IF
     378              :             END IF
     379              :          END DO
     380              : 
     381              :       END DO
     382              : !$OMP END PARALLEL
     383            0 :       CALL o3c_iterator_release(o3c_iterator)
     384              : 
     385            0 :       CALL timestop(handle)
     386              : 
     387            0 :    END SUBROUTINE contract12_o3c
     388              : 
     389              : ! **************************************************************************************************
     390              : !> \brief Contraction of 3-tensor over index 3
     391              : !>        h(ij) = h(ij) + sum_k (ijk)*v(k)
     392              : !> \param o3c ...
     393              : !> \param vec ...
     394              : !> \param matrix ...
     395              : ! **************************************************************************************************
     396            0 :    SUBROUTINE contract3_o3c(o3c, vec, matrix)
     397              :       TYPE(o3c_container_type), POINTER                  :: o3c
     398              :       TYPE(o3c_vec_type), DIMENSION(:), POINTER          :: vec
     399              :       TYPE(dbcsr_type)                                   :: matrix
     400              : 
     401              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract3_o3c'
     402              : 
     403              :       INTEGER                                            :: handle, iatom, icol, ik, irow, jatom, &
     404              :                                                             katom, mepos, nk, nthread, s1, s2
     405              :       LOGICAL                                            :: found, ijsymmetric, trans
     406            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     407            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: v
     408            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pblock
     409            0 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: iabc
     410              :       TYPE(o3c_iterator_type)                            :: o3c_iterator
     411              : 
     412            0 :       CALL timeset(routineN, handle)
     413              : 
     414            0 :       CALL get_o3c_container(o3c, ijsymmetric=ijsymmetric)
     415            0 :       CPASSERT(ijsymmetric)
     416              : 
     417              :       nthread = 1
     418            0 : !$    nthread = omp_get_max_threads()
     419            0 :       CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
     420              : 
     421              : !$OMP PARALLEL DEFAULT(NONE) &
     422              : !$OMP SHARED (nthread,o3c_iterator,vec,matrix)&
     423            0 : !$OMP PRIVATE (mepos,iabc,iatom,jatom,katom,irow,icol,trans,pblock,v,found,ik,nk,work,s1,s2)
     424              : 
     425              :       mepos = 0
     426              : !$    mepos = omp_get_thread_num()
     427              : 
     428              :       DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
     429              :          CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, iatom=iatom, jatom=jatom, katom=katom, &
     430              :                                     integral=iabc)
     431              : 
     432              :          CALL get_o3c_vec(vec, katom, v)
     433              :          nk = SIZE(v)
     434              : 
     435              :          IF (iatom <= jatom) THEN
     436              :             irow = iatom
     437              :             icol = jatom
     438              :             trans = .FALSE.
     439              :          ELSE
     440              :             irow = jatom
     441              :             icol = iatom
     442              :             trans = .TRUE.
     443              :          END IF
     444              : 
     445              :          CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, BLOCK=pblock, found=found)
     446              : 
     447              :          IF (found) THEN
     448              :             s1 = SIZE(pblock, 1); s2 = SIZE(pblock, 2)
     449              :             ALLOCATE (work(s1, s2))
     450              :             work(:, :) = 0.0_dp
     451              : 
     452              :             IF (trans) THEN
     453              :                DO ik = 1, nk
     454              :                   CALL daxpy(s1*s2, v(ik), TRANSPOSE(iabc(:, :, ik)), 1, work(:, :), 1)
     455              :                END DO
     456              :             ELSE
     457              :                DO ik = 1, nk
     458              :                   CALL daxpy(s1*s2, v(ik), iabc(:, :, ik), 1, work(:, :), 1)
     459              :                END DO
     460              :             END IF
     461              : 
     462              :             ! Multiple threads with same irow, icol but different katom (same even in PBCs) can try
     463              :             ! to access the dbcsr block at the same time. Prevent that by CRITICAL section but keep
     464              :             ! computations before hand in order to retain speed
     465              : 
     466              : !$OMP CRITICAL
     467              :             CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, BLOCK=pblock, found=found)
     468              :             CALL daxpy(s1*s2, 1.0_dp, work(:, :), 1, pblock(:, :), 1)
     469              : !$OMP END CRITICAL
     470              : 
     471              :             DEALLOCATE (work)
     472              :          END IF
     473              : 
     474              :       END DO
     475              : !$OMP END PARALLEL
     476            0 :       CALL o3c_iterator_release(o3c_iterator)
     477              : 
     478            0 :       CALL timestop(handle)
     479              : 
     480            0 :    END SUBROUTINE contract3_o3c
     481              : 
     482              : END MODULE qs_o3c_methods
        

Generated by: LCOV version 2.0-1