Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculation of the angular momentum integrals over
      10              : !>      Cartesian Gaussian-type functions.
      11              : !> \par Literature
      12              : !>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
      13              : !> \par History
      14              : !>      none
      15              : !> \author JGH (20.02.2005)
      16              : ! **************************************************************************************************
      17              : MODULE ai_angmom
      18              : 
      19              :    USE ai_cartesian_operators,          ONLY: apply_cartesian_1d
      20              :    USE ai_os_rr,                        ONLY: os_rr_ovlp
      21              :    USE kinds,                           ONLY: dp
      22              :    USE mathconstants,                   ONLY: pi
      23              :    USE orbital_pointers,                ONLY: indco,&
      24              :                                               ncoset
      25              : #include "../base/base_uses.f90"
      26              : 
      27              :    IMPLICIT NONE
      28              : 
      29              :    PRIVATE
      30              : 
      31              : ! *** Public subroutines ***
      32              : 
      33              :    PUBLIC :: angmom
      34              : 
      35              : CONTAINS
      36              : 
      37              : ! **************************************************************************************************
      38              : !> \brief Real angular-momentum kernel (r-C) cross grad, equal to i times L.
      39              : !> \param la_max ...
      40              : !> \param npgfa ...
      41              : !> \param zeta ...
      42              : !> \param rpgfa ...
      43              : !> \param la_min lowest occupied shell; the output retains the full cumulative Cartesian layout
      44              : !> \param lb_max ...
      45              : !> \param npgfb ...
      46              : !> \param zetb ...
      47              : !> \param rpgfb ...
      48              : !> \param rac ...
      49              : !> \param rbc ...
      50              : !> \param angab ...
      51              : ! **************************************************************************************************
      52      1040345 :    SUBROUTINE angmom(la_max, npgfa, zeta, rpgfa, la_min, &
      53       208069 :                      lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
      54              :       INTEGER, INTENT(IN)                                :: la_max, npgfa
      55              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta, rpgfa
      56              :       INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
      57              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zetb, rpgfb
      58              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac, rbc
      59              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: angab
      60              : 
      61              :       INTEGER                                            :: avec(3), bvec(3), coa, cob, ia, ib, &
      62              :                                                             idir, ipgf, jdir, jpgf, kdir, ldrr, &
      63              :                                                             ma, mb
      64              :       REAL(KIND=dp)                                      :: a, b, dab, f0, rab(3), rab2, rap(3), &
      65              :                                                             rbp(3), zet
      66              :       REAL(KIND=dp), ALLOCATABLE                         :: drr(:, :, :), position(:, :, :), &
      67              :                                                             rr(:, :, :)
      68              : 
      69              :       MARK_USED(la_min)
      70       832276 :       rab = rbc - rac
      71       832276 :       dab = NORM2(rab)
      72       208069 :       rab2 = dab**2
      73       208069 :       ldrr = MAX(la_max, lb_max) + 2
      74              :       ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), position(0:ldrr - 1, 0:ldrr - 1, 3), &
      75      1872621 :                 drr(0:ldrr - 1, 0:ldrr - 1, 3))
      76    135937480 :       angab = 0.0_dp
      77       719070 :       DO ipgf = 1, npgfa
      78       511001 :          ma = (ipgf - 1)*ncoset(la_max)
      79      2100482 :          DO jpgf = 1, npgfb
      80      1381412 :             IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
      81       532652 :             mb = (jpgf - 1)*ncoset(lb_max)
      82       532652 :             a = zeta(ipgf)
      83       532652 :             b = zetb(jpgf)
      84       532652 :             zet = a + b
      85      2130608 :             rap = b*rab/zet
      86      2130608 :             rbp = -a*rab/zet
      87       532652 :             f0 = (pi/zet)**1.5_dp*EXP(-a*b*rab2/zet)
      88       532652 :             CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max + 1, zet, ldrr, rr)
      89       532652 :             CALL apply_cartesian_1d(rr, position, la_max, lb_max, .FALSE., .FALSE., 0.0_dp, rac)
      90       532652 :             CALL apply_cartesian_1d(rr, drr, la_max, lb_max, .TRUE., .TRUE., b, [0.0_dp, 0.0_dp, 0.0_dp])
      91              : 
      92              :             ! The ket electronic derivative is the negative of its center derivative.
      93              :             ! Form only the three cross-product components before contraction.
      94      2934400 :             DO cob = 1, ncoset(lb_max)
      95      7562988 :                bvec = indco(:, cob)
      96      1890747 :                ib = mb + cob
      97     11077631 :                DO coa = 1, ncoset(la_max)
      98     31221888 :                   avec = indco(:, coa)
      99      7805472 :                   ia = ma + coa
     100     33112635 :                   DO idir = 1, 3
     101     23416416 :                      jdir = MOD(idir, 3) + 1
     102     23416416 :                      kdir = MOD(idir + 1, 3) + 1
     103              :                      angab(ia, ib, idir) = f0*rr(avec(idir), bvec(idir), idir)* &
     104              :                                            (position(avec(kdir), bvec(kdir), kdir)*drr(avec(jdir), bvec(jdir), jdir) - &
     105     31221888 :                                             position(avec(jdir), bvec(jdir), jdir)*drr(avec(kdir), bvec(kdir), kdir))
     106              :                   END DO
     107              :                END DO
     108              :             END DO
     109              :          END DO
     110              :       END DO
     111       208069 :    END SUBROUTINE angmom
     112              : 
     113              : END MODULE ai_angmom