ai_moments.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Calculation of the moment integrals over Cartesian Gaussian-type
!>      functions.
!> \par Literature
!>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
!> \par History
!>      none
!> \author J. Hutter (16.02.2005)
! **************************************************************************************************
MODULE ai_moments
 
! ax,ay,az  : Angular momentum index numbers of orbital a.
! bx,by,bz  : Angular momentum index numbers of orbital b.
! coset     : Cartesian orbital set pointer.
! dab       : Distance between the atomic centers a and b.
! l{a,b}    : Angular momentum quantum number of shell a or b.
! l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
! l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
! rac       : Distance vector between the atomic center a and reference point c.
! rbc       : Distance vector between the atomic center b and reference point c.
! rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
! zet{a,b}  : Exponents of the Gaussian-type functions a or b.
! zetp      : Reciprocal of the sum of the exponents of orbital a and b.
 
   USE ai_cartesian_operators, ONLY: apply_cartesian_1d, &
                                     apply_cartesian_operator, &
                                     apply_cartesian_operator_pair, &
                                     cartesian_ladder_map, &
                                     cartesian_polynomial_1d
   USE ai_os_rr, ONLY: os_rr_ovlp
   USE ai_overlap, ONLY: overlap_ab
   USE kinds, ONLY: dp
   USE mathconstants, ONLY: pi
   USE orbital_pointers, ONLY: indco, &
                               ncoset
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_moments'
 
   PUBLIC :: expikr, moment, diff_momop, contract_expikr, dipole_force
 
CONTAINS
 
! **************************************************************************************************
!> \brief Contract complex primitive integrals into real cosine and sine AO blocks.
!> \param cos_block ...
!> \param sin_block ...
!> \param iatom ...
!> \param ncoa ...
!> \param nsgfa ...
!> \param sgfa ...
!> \param sphi_a ...
!> \param ldsa ...
!> \param jatom ...
!> \param ncob ...
!> \param nsgfb ...
!> \param sgfb ...
!> \param sphi_b ...
!> \param ldsb ...
!> \param expab complex primitive plane-wave integrals
!> \param work ...
!> \param ldwork ...
!> \param ordered retain the ordered atom pair instead of mapping to an upper block
! **************************************************************************************************
   SUBROUTINE contract_expikr(cos_block, sin_block, &
                              iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, &
                              jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, &
                              expab, work, ldwork, ordered)
 
      REAL(dp), DIMENSION(:, :), POINTER                 :: cos_block, sin_block
      INTEGER, INTENT(IN)                                :: iatom, ncoa, nsgfa, sgfa
      REAL(dp), DIMENSION(:, :), INTENT(IN)              :: sphi_a
      INTEGER, INTENT(IN)                                :: ldsa, jatom, ncob, nsgfb, sgfb
      REAL(dp), DIMENSION(:, :), INTENT(IN)              :: sphi_b
      INTEGER, INTENT(IN)                                :: ldsb
      COMPLEX(dp), DIMENSION(:, :), INTENT(IN)           :: expab
      REAL(dp), DIMENSION(:, :)                          :: work
      INTEGER, INTENT(IN)                                :: ldwork
      LOGICAL, INTENT(IN), OPTIONAL                      :: ordered
 
      LOGICAL                                            :: direct_block
      REAL(dp), ALLOCATABLE                              :: primitive(:, :)
 
      direct_block = iatom <= jatom
      IF (PRESENT(ordered)) direct_block = direct_block .OR. ordered
 
      ALLOCATE (primitive(ncoa, ncob))
            primitive(:, :) = real(expab(1:ncoa, 1:ncob), kind=dp)
         CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
                    1.0_dp, primitive(1, 1), ncoa, &
                    sphi_b(1, sgfb), ldsb, &
                    0.0_dp, work(1, 1), ldwork)
 
         IF (direct_block) THEN
            CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
                       1.0_dp, sphi_a(1, sgfa), ldsa, &
                       work(1, 1), ldwork, &
                       1.0_dp, cos_block(sgfa, sgfb), &
                       SIZE(cos_block, 1))
         ELSE
            CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
                       1.0_dp, work(1, 1), ldwork, &
                       sphi_a(1, sgfa), ldsa, &
                       1.0_dp, cos_block(sgfb, sgfa), &
                       SIZE(cos_block, 1))
         END IF
            primitive(:, :) = aimag(expab(1:ncoa, 1:ncob))
         CALL dgemm("N", "N", ncoa, nsgfb, ncob, &
                    1.0_dp, primitive(1, 1), ncoa, &
                    sphi_b(1, sgfb), ldsb, &
                    0.0_dp, work(1, 1), ldwork)
 
         IF (direct_block) THEN
            CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, &
                       1.0_dp, sphi_a(1, sgfa), ldsa, &
                       work(1, 1), ldwork, &
                       1.0_dp, sin_block(sgfa, sgfb), &
                       SIZE(sin_block, 1))
         ELSE
            CALL dgemm("T", "N", nsgfb, nsgfa, ncoa, &
                       1.0_dp, work(1, 1), ldwork, &
                       sphi_a(1, sgfa), ldsa, &
                       1.0_dp, sin_block(sgfb, sgfa), &
                       SIZE(sin_block, 1))
         END IF
 
   END SUBROUTINE contract_expikr
 
! **************************************************************************************************
!> \brief Evaluate primitive exp(i k.r) integrals through the complex overlap recurrence.
!> \param la_max_set highest bra shell
!> \param npgfa number of bra primitives
!> \param zeta bra Gaussian exponents
!> \param rpgfa bra primitive radii
!> \param la_min_set lowest bra shell
!> \param lb_max highest ket shell
!> \param npgfb number of ket primitives
!> \param zetb ket Gaussian exponents
!> \param rpgfb ket primitive radii
!> \param lb_min lowest ket shell
!> \param rac bra-center coordinates relative to the phase origin
!> \param rbc ket-center coordinates relative to the phase origin
!> \param kvec plane-wave vector
!> \param expab primitive integrals in the full Cartesian layout, zero below the requested shells
!> \param dexpab optional bra-center derivatives with the plane-wave operator held fixed
!> \param dexpba optional ket-center derivatives in the same bra/ket layout
!> \param momentab optional coordinate-weighted integrals relative to the phase origin
!>        Derivative and moment outputs cover the full Cartesian range.
! **************************************************************************************************
   SUBROUTINE expikr(la_max_set, npgfa, zeta, rpgfa, la_min_set, &
                     lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, expab, dexpab, dexpba, momentab)
      INTEGER, INTENT(IN)                                :: la_max_set, npgfa
      REAL(kind=dp), INTENT(IN)                          :: zeta(:), rpgfa(:)
      INTEGER, INTENT(IN)                                :: la_min_set, lb_max, npgfb
      REAL(kind=dp), INTENT(IN)                          :: zetb(:), rpgfb(:)
      INTEGER, INTENT(IN)                                :: lb_min
      REAL(kind=dp), INTENT(IN)                          :: rac(3), rbc(3), kvec(3)
      COMPLEX(KIND=dp), INTENT(INOUT)                    :: expab(:, :)
      COMPLEX(KIND=dp), INTENT(INOUT), OPTIONAL, TARGET  :: dexpab(:, :, :), dexpba(:, :, :), &
                                                            momentab(:, :, :)
 
      COMPLEX(KIND=dp)                                   :: prefactor, zrap(3), zrbp(3)
      COMPLEX(KIND=dp), ALLOCATABLE                      :: opr_rr(:, :, :), rr(:, :, :)
      COMPLEX(KIND=dp), POINTER                          :: opr_values(:, :, :)
      INTEGER                                            :: avec(3), bvec(3), coa, cob, ia, ib, &
                                                            idir, ipgf, jdir, jpgf, kdir, ldrr, &
                                                            lma, lmb, ma, mb, na, nb, op
      LOGICAL                                            :: derivative, right
      REAL(kind=dp)                                      :: dab, exponent, f0, f1, f2, k2, kdp, &
                                                            origin(3), rab(3), rab2, zet, zetp
 
      rab = rbc - rac
      dab = norm2(rab)
      rab2 = sum(rab**2)
      k2 = sum(kvec**2)
      lma = la_max_set
      lmb = lb_max
      IF (PRESENT(dexpab) .OR. PRESENT(momentab)) lma = lma + 1
      IF (PRESENT(dexpba)) lmb = lmb + 1
      ldrr = max(lma, lmb) + 1
      ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3))
      IF (PRESENT(dexpab) .OR. PRESENT(dexpba) .OR. PRESENT(momentab)) THEN
         ALLOCATE (opr_rr(0:la_max_set, 0:lb_max, 3))
      END IF
 
      na = ncoset(la_max_set)
      nb = ncoset(lb_max)
      expab(1:na*npgfa, 1:nb*npgfb) = 0.0_dp
      IF (PRESENT(dexpab)) dexpab(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
      IF (PRESENT(dexpba)) dexpba(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
      IF (PRESENT(momentab)) momentab(1:na*npgfa, 1:nb*npgfb, 1:3) = 0.0_dp
 
      DO ipgf = 1, npgfa
         ma = (ipgf - 1)*na
         DO jpgf = 1, npgfb
            IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
            mb = (jpgf - 1)*nb
            zet = zeta(ipgf) + zetb(jpgf)
            zetp = 1.0_dp/zet
            f0 = (pi*zetp)**1.5_dp
            f1 = zetb(jpgf)*zetp
            f2 = 0.5_dp*zetp
            kdp = zetp*dot_product(kvec, zeta(ipgf)*rac + zetb(jpgf)*rbc)
            prefactor = f0*exp(cmplx(-zeta(ipgf)*f1*rab2 - 0.25_dp*k2*zetp, kdp, kind=dp))
 
            ! The normalized 1D tables carry no three-dimensional prefactor.
            zrap = cmplx(f1*rab, f2*kvec, kind=dp)
            zrbp = zrap - rab
            CALL os_rr_ovlp(zrap, lma, zrbp, lmb, zet, ldrr, rr)
            DO cob = ncoset(lb_min - 1) + 1, nb
               bvec = indco(:, cob)
               DO coa = ncoset(la_min_set - 1) + 1, na
                  avec = indco(:, coa)
                  expab(ma + coa, mb + cob) = prefactor*rr(avec(1), bvec(1), 1)* &
                                              rr(avec(2), bvec(2), 2)*rr(avec(3), bvec(3), 3)
               END DO
            END DO
 
            DO op = 1, 3
               SELECT CASE (op)
               CASE (1)
                  IF (.NOT. PRESENT(dexpab)) cycle
                  opr_values => dexpab
                  exponent = zeta(ipgf)
               CASE (2)
                  IF (.NOT. PRESENT(dexpba)) cycle
                  opr_values => dexpba
                  exponent = zetb(jpgf)
               CASE (3)
                  IF (.NOT. PRESENT(momentab)) cycle
                  opr_values => momentab
                  exponent = 0.0_dp
               END SELECT
               right = op == 2
               derivative = op /= 3
               origin = 0.0_dp
               IF (.NOT. derivative) origin = rac
               CALL apply_cartesian_1d(rr, opr_rr, la_max_set, lb_max, right, derivative, exponent, origin)
               DO cob = 1, nb
                  bvec = indco(:, cob)
                  ib = mb + cob
                  DO coa = 1, na
                     avec = indco(:, coa)
                     ia = ma + coa
                     DO idir = 1, 3
                        jdir = mod(idir, 3) + 1
                        kdir = mod(idir + 1, 3) + 1
                        opr_values(ia, ib, idir) = prefactor*opr_rr(avec(idir), bvec(idir), idir)* &
                                                   rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir)
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
   END SUBROUTINE expikr
 
! **************************************************************************************************
!> \brief ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param la_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param lc_max ...
!> \param rac ...
!> \param rbc ...
!> \param mab ...
!> \param overlap optional zeroth moment in primitive-major cumulative Cartesian order;
!>        entries outside the requested row and column ranges are untouched
! **************************************************************************************************
   SUBROUTINE moment(la_max, npgfa, zeta, rpgfa, la_min, &
                     lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab, overlap)
      INTEGER, INTENT(IN)                                :: la_max, npgfa
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zeta, rpgfa
      INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zetb, rpgfb
      INTEGER, INTENT(IN)                                :: lc_max
      REAL(kind=dp), DIMENSION(3), INTENT(IN)            :: rac, rbc
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: mab
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: overlap
 
      INTEGER                                            :: ax, ay, az, bx, by, bz, coa, cob, i, ia, &
                                                            ib, imom, ipgf, jpgf, ldrr, lx, ly, &
                                                            lz, ma, mb, powers(0:lc_max)
      REAL(kind=dp)                                      :: a, b, dab, f0, rab(3), rab2, rap(3), &
                                                            rbp(3), zet
      REAL(kind=dp), ALLOCATABLE                         :: polynomial(:, :, :, :), rr(:, :, :)
 
      mark_used(la_min)
      rab = rbc - rac
      dab = norm2(rab)
      rab2 = dab**2
      ldrr = max(la_max + lc_max, lb_max) + 1
      ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), polynomial(0:la_max, 0:lb_max, 0:lc_max, 3))
      powers = [(i, i=0, lc_max)]
      mab(:npgfa*ncoset(la_max), :npgfb*ncoset(lb_max), :ncoset(lc_max) - 1) = 0.0_dp
      IF (PRESENT(overlap)) overlap(:npgfa*ncoset(la_max), :npgfb*ncoset(lb_max)) = 0.0_dp
      DO ipgf = 1, npgfa
         ma = (ipgf - 1)*ncoset(la_max)
         DO jpgf = 1, npgfb
            IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) cycle
            mb = (jpgf - 1)*ncoset(lb_max)
            a = zeta(ipgf)
            b = zetb(jpgf)
            zet = a + b
            rap = b*rab/zet
            rbp = -a*rab/zet
            f0 = (pi/zet)**1.5_dp*exp(-a*b*rab2/zet)
            CALL os_rr_ovlp(rap, la_max + lc_max, rbp, lb_max, zet, ldrr, rr)
            CALL cartesian_polynomial_1d(rr, polynomial, la_max, lb_max, powers, rac)
 
            DO cob = 1, ncoset(lb_max)
               bx = indco(1, cob)
               by = indco(2, cob)
               bz = indco(3, cob)
               ib = mb + cob
               DO coa = 1, ncoset(la_max)
                  ax = indco(1, coa)
                  ay = indco(2, coa)
                  az = indco(3, coa)
                  ia = ma + coa
                  IF (PRESENT(overlap)) THEN
                     overlap(ia, ib) = f0*polynomial(ax, bx, 0, 1)* &
                                       polynomial(ay, by, 0, 2)*polynomial(az, bz, 0, 3)
                  END IF
                  ! Preserve the coset ordering and omission of the zeroth moment.
                  DO imom = 1, ncoset(lc_max) - 1
                     lx = indco(1, imom + 1)
                     ly = indco(2, imom + 1)
                     lz = indco(3, imom + 1)
                     mab(ia, ib, imom) = f0*polynomial(ax, bx, lx, 1)* &
                                         polynomial(ay, by, ly, 2)*polynomial(az, bz, lz, 3)
                  END DO
               END DO
            END DO
         END DO
      END DO
   END SUBROUTINE moment
 
! **************************************************************************************************
!> \brief This returns the derivative of the moment integrals [a|\mu|b].
!>       By default, it differentiates the primitive on the right:
!>       [a|\mu|d/dR_bi] =  2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
!>       A weighted derivative combines the left and right primitive derivatives
!>       using deltaR for the corresponding atom centers.
!>       order indicates the max order of the moment operator to be calculated
!>       1: dipole
!>       2: quadrupole
!>       ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param la_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param lb_min ...
!> \param order ...
!> \param rac ...
!> \param rbc ...
!> \param difmab ...
!> \param mab_ext optional precomputed moments with zero screened primitive blocks, including raised shells
!> \param deltaR optional weights for the left and right primitive derivatives
!> \param lambda optional atom selector for the factor (iatom == lambda) - (jatom == lambda)
!> \param iatom atom associated with the left basis function
!> \param jatom atom associated with the right basis function
!> \note
! **************************************************************************************************
   SUBROUTINE diff_momop(la_max, npgfa, zeta, rpgfa, la_min, &
                         lb_max, npgfb, zetb, rpgfb, lb_min, &
                         order, rac, rbc, difmab, mab_ext, deltaR, lambda, iatom, jatom)
 
      INTEGER, INTENT(IN)                                :: la_max, npgfa
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zeta, rpgfa
      INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zetb, rpgfb
      INTEGER, INTENT(IN)                                :: lb_min, order
      REAL(kind=dp), DIMENSION(3), INTENT(IN)            :: rac, rbc
      REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(OUT)  :: difmab
      REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
         POINTER                                         :: mab_ext
      REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: deltar
      INTEGER, INTENT(IN), OPTIONAL                      :: lambda, iatom, jatom
 
      INTEGER                                            :: ider, imom, lda, lda_min, ldb, ldb_min
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: difmab_tmp
      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: mab
 
      lda_min = max(0, la_min - 1)
      ldb_min = max(0, lb_min - 1)
      lda = ncoset(la_max)*npgfa
      ldb = ncoset(lb_max)*npgfb
      ALLOCATE (difmab_tmp(lda, ldb, 3))
 
      IF (PRESENT(mab_ext)) THEN
         mab => mab_ext
      ELSE
         ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
                       ncoset(order) - 1))
         mab = 0.0_dp
!     *** Calculate the primitive overlap integrals ***
         CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
                     lb_max + 1, npgfb, zetb, rpgfb, &
                     order, rac, rbc, mab)
 
      END IF
      DO imom = 1, ncoset(order) - 1
         CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, mab(:, :, imom), &
                                       difmab_tmp, right=.true., derivative=.true., exponents=zetb)
         IF (PRESENT(deltar)) THEN
            cpassert(ASSOCIATED(deltar))
            cpassert(PRESENT(iatom) .AND. PRESENT(jatom))
            DO ider = 1, 3
               difmab(1:lda, 1:ldb, imom, ider) = &
                  difmab_tmp(1:lda, 1:ldb, ider)*deltar(ider, jatom)
            END DO
 
            CALL apply_cartesian_operator(la_max, 0, npgfa, lb_max, 0, npgfb, mab(:, :, imom), &
                                          difmab_tmp, right=.false., derivative=.true., exponents=zeta)
            DO ider = 1, 3
               difmab(1:lda, 1:ldb, imom, ider) = difmab(1:lda, 1:ldb, imom, ider) &
                                                  + difmab_tmp(1:lda, 1:ldb, ider)*deltar(ider, iatom)
            END DO
         ELSE
            difmab(1:lda, 1:ldb, imom, :) = difmab_tmp(1:lda, 1:ldb, :)
         END IF
      END DO
 
      IF (PRESENT(lambda)) THEN
         cpassert(.NOT. PRESENT(deltar))
         cpassert(PRESENT(iatom) .AND. PRESENT(jatom))
         IF (iatom == lambda .AND. jatom == lambda) THEN
            difmab = 0.0_dp
         ELSE IF (iatom == lambda) THEN
            ! The right-hand derivative is selected with a positive sign.
         ELSE IF (jatom == lambda) THEN
            difmab = -difmab
         ELSE
            difmab = 0.0_dp
         END IF
      END IF
 
      IF (PRESENT(mab_ext)) THEN
         NULLIFY (mab)
      ELSE
         DEALLOCATE (mab)
      END IF
      DEALLOCATE (difmab_tmp)
 
   END SUBROUTINE diff_momop
 
! **************************************************************************************************
!> \brief This returns the derivative of the dipole integrals [a|x|b], with respect
!>       to the position of the primitive on the left and right, i.e.
!>       [da/dR_ai|\mu|b] =  2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param la_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param lb_min ...
!> \param order ...
!> \param rac ...
!> \param rbc ...
!> \param pab ...
!> \param forcea ...
!> \param forceb ...
!> \note
! **************************************************************************************************
   SUBROUTINE dipole_force(la_max, npgfa, zeta, rpgfa, la_min, &
                           lb_max, npgfb, zetb, rpgfb, lb_min, &
                           order, rac, rbc, pab, forcea, forceb)
 
      INTEGER, INTENT(IN)                                :: la_max, npgfa
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zeta, rpgfa
      INTEGER, INTENT(IN)                                :: la_min, lb_max, npgfb
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: zetb, rpgfb
      INTEGER, INTENT(IN)                                :: lb_min, order
      REAL(kind=dp), DIMENSION(3), INTENT(IN)            :: rac, rbc
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: pab
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), TARGET :: forcea, forceb
 
      INTEGER, PARAMETER :: components(2, 9) = reshape([1, 1, 2, 1, 3, 1, 1, 2, 2, 2, 3, 2, 1, 3, 2, 3, 3, 3], [2, 9])
      INTEGER :: axis, first_a, first_b, icomp, ipgf, jpgf, ma, mb, nrow, ncol, side
      INTEGER :: map_a(4, ncoset(la_max) - ncoset(la_min - 1), 3), &
                 map_b(4, ncoset(lb_max) - ncoset(lb_min - 1), 3)
      REAL(kind=dp) :: origin(3, 2), rab(3), distance, rr_work(3*(max(la_max, lb_max) + 2)**2)
      REAL(kind=dp), ALLOCATABLE :: source(:, :), work(:, :), derivatives(:, :, :)
      REAL(kind=dp), POINTER :: forces(:, :)
 
      cpassert(order == 1)
      mark_used(order)
      first_a = ncoset(la_min - 1) + 1
      first_b = ncoset(lb_min - 1) + 1
      nrow = ncoset(la_max) - first_a + 1
      ncol = ncoset(lb_max) - first_b + 1
      ALLOCATE (source(ncoset(la_max + 1), ncoset(lb_max + 1)), work(ncoset(la_max + 1), ncoset(lb_max + 1)), &
                derivatives(nrow, ncol, 9))
      DO axis = 1, 3
         CALL cartesian_ladder_map(la_max, la_min, axis, map_a(:, :, axis), source_lmin=0)
         CALL cartesian_ladder_map(lb_max, lb_min, axis, map_b(:, :, axis), source_lmin=0)
      END DO
      origin = 0.0_dp
      rab = rbc - rac
      distance = norm2(rab)
      DO ipgf = 1, npgfa
         ma = (ipgf - 1)*ncoset(la_max)
         DO jpgf = 1, npgfb
            IF (rpgfa(ipgf) + rpgfb(jpgf) < distance) cycle
            mb = (jpgf - 1)*ncoset(lb_max)
            CALL overlap_ab(la_max + 1, 0, 1, rpgfa(ipgf:ipgf), zeta(ipgf:ipgf), &
                            lb_max + 1, 0, 1, rpgfb(jpgf:jpgf), zetb(jpgf:jpgf), rab, &
                            sab=source, rr_work=rr_work)
            ! Multiply by the coordinate on the side opposite the differentiated Gaussian.
            DO side = 1, 2
               IF (side == 1) THEN
                  forces => forcea
                  origin(:, 1) = rbc
                  CALL apply_cartesian_operator_pair(source, derivatives, map_b, map_a, components, &
                                                     [.false., .true.], [0.0_dp, zeta(ipgf)], &
                                                     [.true., .false.], work(:, 1:ncol), origin)
               ELSE
                  forces => forceb
                  origin(:, 1) = rac
                  CALL apply_cartesian_operator_pair(source, derivatives, map_a, map_b, components, &
                                                     [.false., .true.], [0.0_dp, zetb(jpgf)], &
                                                     [.false., .true.], work(1:nrow, :), origin)
               END IF
               DO icomp = 1, 9
                  forces(components(1, icomp), components(2, icomp)) = &
                     forces(components(1, icomp), components(2, icomp)) + &
                     sum(pab(ma + first_a:ma + ncoset(la_max), mb + first_b:mb + ncoset(lb_max))*derivatives(:, :, icomp))
               END DO
            END DO
         END DO
      END DO
 
   END SUBROUTINE dipole_force
 
END MODULE ai_moments