LCOV - code coverage report
Current view: top level - src - commutator_rkinetic.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 63 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 2 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              : ! **************************************************************************************************
       9              : !> \brief Calculation of commutator of kinetic energy and position operator
      10              : !> \par History
      11              : !>      JGH: from qs_kinetic
      12              : !> \author Juerg Hutter
      13              : ! **************************************************************************************************
      14              : MODULE commutator_rkinetic
      15              :    USE ai_contraction,                  ONLY: block_add,&
      16              :                                               contraction
      17              :    USE ai_kinetic,                      ONLY: kinetic
      18              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      19              :                                               gto_basis_set_type
      20              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
      21              :                                               dbcsr_p_type
      22              :    USE kinds,                           ONLY: dp
      23              :    USE orbital_pointers,                ONLY: coset,&
      24              :                                               ncoset
      25              :    USE qs_integral_utils,               ONLY: basis_set_list_setup,&
      26              :                                               get_memory_usage
      27              :    USE qs_kind_types,                   ONLY: qs_kind_type
      28              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      29              :                                               get_neighbor_list_set_p,&
      30              :                                               neighbor_list_iterate,&
      31              :                                               neighbor_list_iterator_create,&
      32              :                                               neighbor_list_iterator_p_type,&
      33              :                                               neighbor_list_iterator_release,&
      34              :                                               neighbor_list_set_p_type
      35              : 
      36              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
      37              : #include "./base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              : 
      41              :    PRIVATE
      42              : 
      43              : ! *** Global parameters ***
      44              : 
      45              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rkinetic'
      46              : 
      47              : ! *** Public subroutines ***
      48              : 
      49              :    PUBLIC :: build_com_tr_matrix
      50              : 
      51              : CONTAINS
      52              : 
      53              : ! **************************************************************************************************
      54              : !> \brief   Calculation of commutator [T,r] over Cartesian Gaussian functions.
      55              : !> \param   matrix_tr ...
      56              : !> \param   qs_kind_set ...
      57              : !> \param   basis_type basis set to be used
      58              : !> \param   sab_nl pair list (must be consistent with basis sets!)
      59              : !> \date    11.10.2010
      60              : !> \par     History
      61              : !>          Ported from qs_overlap, replaces code in build_core_hamiltonian
      62              : !>          Refactoring [07.2014] JGH
      63              : !>          Simplify options and use new kinetic energy integral routine
      64              : !>          Adapted from qs_kinetic [07.2016]
      65              : !> \author  JGH
      66              : !> \version 1.0
      67              : ! **************************************************************************************************
      68            0 :    SUBROUTINE build_com_tr_matrix(matrix_tr, qs_kind_set, basis_type, sab_nl)
      69              : 
      70              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_tr
      71              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      72              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
      73              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      74              :          POINTER                                         :: sab_nl
      75              : 
      76              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_com_tr_matrix'
      77              : 
      78              :       INTEGER                                            :: handle, iatom, icol, ikind, ir, irow, &
      79              :                                                             iset, jatom, jkind, jset, ldsab, ltab, &
      80              :                                                             mepos, ncoa, ncob, nkind, nseta, &
      81              :                                                             nsetb, nthread, sgfa, sgfb
      82              :       INTEGER, DIMENSION(3)                              :: cell
      83            0 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
      84            0 :                                                             npgfb, nsgfa, nsgfb
      85            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
      86              :       LOGICAL                                            :: do_symmetric, found, trans
      87              :       REAL(KIND=dp)                                      :: rab2, tab
      88            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: qab, tkab
      89            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: kab
      90              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
      91            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
      92            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: kx_block, ky_block, kz_block, rpgfa, &
      93            0 :                                                             rpgfb, scon_a, scon_b, zeta, zetb
      94            0 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
      95              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
      96              :       TYPE(neighbor_list_iterator_p_type), &
      97            0 :          DIMENSION(:), POINTER                           :: nl_iterator
      98              : 
      99            0 :       CALL timeset(routineN, handle)
     100              : 
     101            0 :       nkind = SIZE(qs_kind_set)
     102              : 
     103              :       ! check for symmetry
     104            0 :       CPASSERT(SIZE(sab_nl) > 0)
     105            0 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
     106              : 
     107              :       ! prepare basis set
     108            0 :       ALLOCATE (basis_set_list(nkind))
     109            0 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
     110              : 
     111              :       ! *** Allocate work storage ***
     112            0 :       ldsab = get_memory_usage(qs_kind_set, basis_type)
     113              : 
     114              :       nthread = 1
     115            0 : !$    nthread = omp_get_max_threads()
     116              :       ! Iterate of neighbor list
     117            0 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
     118              : 
     119              : !$OMP PARALLEL DEFAULT(NONE) &
     120              : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric,&
     121              : !$OMP         ncoset,matrix_tr,basis_set_list)  &
     122              : !$OMP PRIVATE (kx_block,ky_block,kz_block,mepos,kab,qab,tab,ikind,jkind,iatom,jatom,rab,cell,&
     123              : !$OMP          basis_set_a,basis_set_b,&
     124              : !$OMP          first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, &
     125              : !$OMP          zeta, first_sgfb, lb_max, lb_min, ltab, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, tkab, &
     126            0 : !$OMP          zetb, scon_a, scon_b, irow, icol, found, trans, rab2, sgfa, sgfb, iset, jset)
     127              : 
     128              :       mepos = 0
     129              : !$    mepos = omp_get_thread_num()
     130              : 
     131              :       ALLOCATE (kab(ldsab, ldsab, 3), qab(ldsab, ldsab))
     132              : 
     133              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     134              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
     135              :                                 iatom=iatom, jatom=jatom, r=rab, cell=cell)
     136              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     137              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     138              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     139              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     140              :          ! basis ikind
     141              :          first_sgfa => basis_set_a%first_sgf
     142              :          la_max => basis_set_a%lmax
     143              :          la_min => basis_set_a%lmin
     144              :          npgfa => basis_set_a%npgf
     145              :          nseta = basis_set_a%nset
     146              :          nsgfa => basis_set_a%nsgf_set
     147              :          rpgfa => basis_set_a%pgf_radius
     148              :          set_radius_a => basis_set_a%set_radius
     149              :          scon_a => basis_set_a%scon
     150              :          zeta => basis_set_a%zet
     151              :          ! basis jkind
     152              :          first_sgfb => basis_set_b%first_sgf
     153              :          lb_max => basis_set_b%lmax
     154              :          lb_min => basis_set_b%lmin
     155              :          npgfb => basis_set_b%npgf
     156              :          nsetb = basis_set_b%nset
     157              :          nsgfb => basis_set_b%nsgf_set
     158              :          rpgfb => basis_set_b%pgf_radius
     159              :          set_radius_b => basis_set_b%set_radius
     160              :          scon_b => basis_set_b%scon
     161              :          zetb => basis_set_b%zet
     162              : 
     163              :          IF (do_symmetric) THEN
     164              :             IF (iatom <= jatom) THEN
     165              :                irow = iatom
     166              :                icol = jatom
     167              :             ELSE
     168              :                irow = jatom
     169              :                icol = iatom
     170              :             END IF
     171              :          ELSE
     172              :             irow = iatom
     173              :             icol = jatom
     174              :          END IF
     175              :          NULLIFY (kx_block)
     176              :          CALL dbcsr_get_block_p(matrix=matrix_tr(1)%matrix, &
     177              :                                 row=irow, col=icol, BLOCK=kx_block, found=found)
     178              :          CPASSERT(found)
     179              :          NULLIFY (ky_block)
     180              :          CALL dbcsr_get_block_p(matrix=matrix_tr(2)%matrix, &
     181              :                                 row=irow, col=icol, BLOCK=ky_block, found=found)
     182              :          CPASSERT(found)
     183              :          NULLIFY (kz_block)
     184              :          CALL dbcsr_get_block_p(matrix=matrix_tr(3)%matrix, &
     185              :                                 row=irow, col=icol, BLOCK=kz_block, found=found)
     186              :          CPASSERT(found)
     187              : 
     188              :          rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     189              :          tab = SQRT(rab2)
     190              :          trans = do_symmetric .AND. (iatom > jatom)
     191              : 
     192              :          DO iset = 1, nseta
     193              : 
     194              :             ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     195              :             sgfa = first_sgfa(1, iset)
     196              : 
     197              :             DO jset = 1, nsetb
     198              : 
     199              :                IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
     200              : 
     201              :                ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     202              :                sgfb = first_sgfb(1, jset)
     203              : 
     204              :                ! calclulate integrals
     205              :                ltab = MAX(npgfa(iset)*ncoset(la_max(iset) + 1), npgfb(jset)*ncoset(lb_max(jset) + 1))
     206              :                ALLOCATE (tkab(ltab, ltab))
     207              :                CALL kinetic(la_max(iset) + 1, la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     208              :                             lb_max(jset) + 1, lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     209              :                             rab, tkab)
     210              :                ! commutator
     211              :                CALL comab_opr(la_max(iset), npgfa(iset), rpgfa(:, iset), la_min(iset), &
     212              :                               lb_max(jset), npgfb(jset), rpgfb(:, jset), lb_min(jset), &
     213              :                               tab, tkab, kab)
     214              :                DEALLOCATE (tkab)
     215              :                ! Contraction step
     216              :                DO ir = 1, 3
     217              :                   CALL contraction(kab(:, :, ir), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
     218              :                                    cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), trans=trans)
     219              : !$OMP CRITICAL(blockadd)
     220              :                   SELECT CASE (ir)
     221              :                   CASE (1)
     222              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kx_block, sgfa, sgfb, trans=trans)
     223              :                   CASE (2)
     224              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), ky_block, sgfa, sgfb, trans=trans)
     225              :                   CASE (3)
     226              :                      CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), kz_block, sgfa, sgfb, trans=trans)
     227              :                   END SELECT
     228              : !$OMP END CRITICAL(blockadd)
     229              :                END DO
     230              : 
     231              :             END DO
     232              :          END DO
     233              : 
     234              :       END DO
     235              :       DEALLOCATE (kab, qab)
     236              : !$OMP END PARALLEL
     237            0 :       CALL neighbor_list_iterator_release(nl_iterator)
     238              : 
     239              :       ! Release work storage
     240            0 :       DEALLOCATE (basis_set_list)
     241              : 
     242            0 :       CALL timestop(handle)
     243              : 
     244            0 :    END SUBROUTINE build_com_tr_matrix
     245              : 
     246              : ! **************************************************************************************************
     247              : !> \brief   Calculate the commutator [O,r] from the integrals [a|O|b].
     248              : !>          We assume that on input all integrals [a+1|O|b+1] are available.
     249              : !>          [a|[O,ri]|b] = [a|O|b+1i] - [a+1i|O|b]
     250              : !> \param la_max ...
     251              : !> \param npgfa ...
     252              : !> \param rpgfa ...
     253              : !> \param la_min ...
     254              : !> \param lb_max ...
     255              : !> \param npgfb ...
     256              : !> \param rpgfb ...
     257              : !> \param lb_min ...
     258              : !> \param dab ...
     259              : !> \param ab ...
     260              : !> \param comabr ...
     261              : !>
     262              : !> \date    25.07.2016
     263              : !> \par Literature
     264              : !>          S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
     265              : !> \par Parameters
     266              : !>      - ax,ay,az  : Angular momentum index numbers of orbital a.
     267              : !>      - bx,by,bz  : Angular momentum index numbers of orbital b.
     268              : !>      - coset     : Cartesian orbital set pointer.
     269              : !>      - l{a,b}    : Angular momentum quantum number of shell a or b.
     270              : !>      - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
     271              : !>      - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
     272              : !>      - ncoset    : Number of orbitals in a Cartesian orbital set.
     273              : !>      - npgf{a,b} : Degree of contraction of shell a or b.
     274              : !>      - rab       : Distance vector between the atomic centers a and b.
     275              : !>      - rab2      : Square of the distance between the atomic centers a and b.
     276              : !>      - rac       : Distance vector between the atomic centers a and c.
     277              : !>      - rac2      : Square of the distance between the atomic centers a and c.
     278              : !>      - rbc       : Distance vector between the atomic centers b and c.
     279              : !>      - rbc2      : Square of the distance between the atomic centers b and c.
     280              : !>      - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
     281              : !>      - zet{a,b}  : Exponents of the Gaussian-type functions a or b.
     282              : !>      - zetp      : Reciprocal of the sum of the exponents of orbital a and b.
     283              : !>
     284              : !> \author  JGH
     285              : !> \version 1.0
     286              : ! **************************************************************************************************
     287            0 :    SUBROUTINE comab_opr(la_max, npgfa, rpgfa, la_min, lb_max, npgfb, rpgfb, lb_min, &
     288            0 :                         dab, ab, comabr)
     289              :       INTEGER, INTENT(IN)                                :: la_max, npgfa
     290              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfa
     291              :       INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
     292              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rpgfb
     293              :       INTEGER, INTENT(IN)                                :: lb_min
     294              :       REAL(KIND=dp), INTENT(IN)                          :: dab
     295              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: ab
     296              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: comabr
     297              : 
     298              :       INTEGER                                            :: ax, ay, az, bx, by, bz, coa, coap, &
     299              :                                                             coapx, coapy, coapz, cob, cobp, cobpx, &
     300              :                                                             cobpy, cobpz, ipgf, jpgf, la, lb, na, &
     301              :                                                             nap, nb, nbp, ofa, ofb
     302              : 
     303            0 :       comabr = 0.0_dp
     304              : 
     305            0 :       ofa = ncoset(la_min - 1)
     306            0 :       ofb = ncoset(lb_min - 1)
     307              : 
     308            0 :       na = 0
     309            0 :       nap = 0
     310            0 :       DO ipgf = 1, npgfa
     311            0 :          nb = 0
     312            0 :          nbp = 0
     313            0 :          DO jpgf = 1, npgfb
     314            0 :             IF (rpgfa(ipgf) + rpgfb(jpgf) > dab) THEN
     315            0 :                DO la = la_min, la_max
     316            0 :                   DO ax = 0, la
     317            0 :                      DO ay = 0, la - ax
     318            0 :                         az = la - ax - ay
     319            0 :                         coa = na + coset(ax, ay, az) - ofa
     320            0 :                         coap = nap + coset(ax, ay, az) - ofa
     321            0 :                         coapx = nap + coset(ax + 1, ay, az) - ofa
     322            0 :                         coapy = nap + coset(ax, ay + 1, az) - ofa
     323            0 :                         coapz = nap + coset(ax, ay, az + 1) - ofa
     324            0 :                         DO lb = lb_min, lb_max
     325            0 :                            DO bx = 0, lb
     326            0 :                               DO by = 0, lb - bx
     327            0 :                                  bz = lb - bx - by
     328            0 :                                  cob = nb + coset(bx, by, bz) - ofb
     329            0 :                                  cobp = nbp + coset(bx, by, bz) - ofb
     330            0 :                                  cobpx = nbp + coset(bx + 1, by, bz) - ofb
     331            0 :                                  cobpy = nbp + coset(bx, by + 1, bz) - ofb
     332            0 :                                  cobpz = nbp + coset(bx, by, bz + 1) - ofb
     333              :                                  ! [a|[O,ri]|b] = [a|O|b+1i] - [a+1i|O|b]
     334            0 :                                  comabr(coa, cob, 1) = ab(coap, cobpx) - ab(coapx, cobp)
     335            0 :                                  comabr(coa, cob, 2) = ab(coap, cobpy) - ab(coapy, cobp)
     336            0 :                                  comabr(coa, cob, 3) = ab(coap, cobpz) - ab(coapz, cobp)
     337              :                               END DO
     338              :                            END DO
     339              :                         END DO
     340              :                      END DO
     341              :                   END DO
     342              :                END DO
     343              :             END IF
     344            0 :             nb = nb + ncoset(lb_max) - ofb
     345            0 :             nbp = nbp + ncoset(lb_max + 1) - ofb
     346              :          END DO
     347            0 :          na = na + ncoset(la_max) - ofa
     348            0 :          nap = nap + ncoset(la_max + 1) - ofa
     349              :       END DO
     350              : 
     351            0 :    END SUBROUTINE comab_opr
     352              : 
     353              : END MODULE commutator_rkinetic
     354              : 
        

Generated by: LCOV version 2.0-1