LCOV - code coverage report
Current view: top level - src - qs_o3c_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 0 53 0.0 %
Date: 2024-03-29 07:50:05 Functions: 0 3 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : ! **************************************************************************************************
       8             : !> \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 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 1.15