pao_linpot_rotinv.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 Rotationally invariant parametrization of Fock matrix.
!> \author Ole Schuett
! **************************************************************************************************
MODULE pao_linpot_rotinv
   USE ai_overlap,                      ONLY: overlap_aab
   USE atomic_kind_types,               ONLY: get_atomic_kind
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: gamma1
   USE mathlib,                         ONLY: multinomial
   USE orbital_pointers,                ONLY: indco,&
                                              ncoset
   USE particle_types,                  ONLY: particle_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              pao_potential_type,&
                                              qs_kind_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_linpot_rotinv'
 
   PUBLIC :: linpot_rotinv_count_terms, linpot_rotinv_calc_terms, linpot_rotinv_calc_forces
 
CONTAINS
 
! **************************************************************************************************
!> \brief Count number of terms for given atomic kind
!> \param qs_env ...
!> \param ikind ...
!> \param nterms ...
! **************************************************************************************************
   SUBROUTINE linpot_rotinv_count_terms(qs_env, ikind, nterms)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      INTEGER, INTENT(IN)                                :: ikind
      INTEGER, INTENT(OUT)                               :: nterms
 
      CHARACTER(len=*), PARAMETER :: routinen = 'linpot_rotinv_count_terms'
 
      INTEGER                                            :: handle, ipot, iset, ishell, ishell_abs, &
                                                            lmax, lmin, lpot, max_shell, &
                                                            min_shell, npots, nshells, pot_maxl
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: shell_l
      TYPE(gto_basis_set_type), POINTER                  :: basis_set
      TYPE(pao_potential_type), DIMENSION(:), POINTER    :: pao_potentials
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
      CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=pao_potentials)
 
      nshells = sum(basis_set%nshell)
      npots = SIZE(pao_potentials)
 
      cpwarn_if(npots == 0, "Found no PAO_POTENTIAL section")
 
      ! fill shell_l
      ALLOCATE (shell_l(nshells))
      DO iset = 1, basis_set%nset
      DO ishell = 1, basis_set%nshell(iset)
         ishell_abs = sum(basis_set%nshell(1:iset - 1)) + ishell
         shell_l(ishell_abs) = basis_set%l(ishell, iset)
      END DO
      END DO
 
      nterms = 0
 
      ! terms sensing neighboring atoms
      DO ipot = 1, npots
         pot_maxl = pao_potentials(ipot)%maxl ! maxl is taken from central atom
         IF (pot_maxl < 0) THEN
            cpabort("ROTINV parametrization requires non-negative PAO_POTENTIAL%MAXL")
         END IF
         IF (mod(pot_maxl, 2) /= 0) THEN
            cpabort("ROTINV parametrization requires even-numbered PAO_POTENTIAL%MAXL")
         END IF
         DO max_shell = 1, nshells
         DO min_shell = 1, max_shell
         DO lpot = 0, pot_maxl, 2
            lmin = shell_l(min_shell)
            lmax = shell_l(max_shell)
            IF (lmin == 0 .AND. lmax == 0) cycle ! covered by central terms
            nterms = nterms + 1
         END DO
         END DO
         END DO
      END DO
 
      ! spherical symmetric terms on central atom
      DO max_shell = 1, nshells
      DO min_shell = 1, max_shell
         IF (shell_l(min_shell) /= shell_l(max_shell)) cycle ! need quadratic block
         nterms = nterms + 1
      END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE linpot_rotinv_count_terms
 
! **************************************************************************************************
!> \brief Calculate all potential terms of the rotinv parametrization
!> \param qs_env ...
!> \param iatom ...
!> \param V_blocks ...
! **************************************************************************************************
   SUBROUTINE linpot_rotinv_calc_terms(qs_env, iatom, V_blocks)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      INTEGER, INTENT(IN)                                :: iatom
      REAL(dp), DIMENSION(:, :, :), INTENT(OUT), TARGET  :: v_blocks
 
      CHARACTER(len=*), PARAMETER :: routinen = 'linpot_rotinv_calc_terms'
 
      INTEGER :: handle, i, ic, ikind, ipot, iset, ishell, ishell_abs, jatom, jkind, jset, jshell, &
         jshell_abs, kterm, la1_max, la1_min, la2_max, la2_min, lb_max, lb_min, lpot, n, na1, na2, &
         natoms, nb, ncfga1, ncfga2, ncfgb, npgfa1, npgfa2, npgfb, npots, pot_maxl, sgfa1, sgfa2, &
         sgla1, sgla2
      REAL(dp)                                           :: coeff, norm2, pot_beta, pot_weight, &
                                                            rpgfa_max, tab
      REAL(dp), DIMENSION(3)                             :: ra, rab, rb
      REAL(dp), DIMENSION(:), POINTER                    :: rpgfa1, rpgfa2, rpgfb, zeta1, zeta2, zetb
      REAL(dp), DIMENSION(:, :), POINTER                 :: t1, t2, v12, v21
      REAL(dp), DIMENSION(:, :, :), POINTER              :: block_v_full, saab, saal
      TYPE(cell_type), POINTER                           :: cell
      TYPE(gto_basis_set_type), POINTER                  :: basis_set
      TYPE(pao_potential_type), DIMENSION(:), POINTER    :: ipao_potentials, jpao_potentials
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, &
                      natom=natoms, &
                      cell=cell, &
                      particle_set=particle_set, &
                      qs_kind_set=qs_kind_set)
 
      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
      CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=ipao_potentials)
      npots = SIZE(ipao_potentials)
      n = basis_set%nsgf ! primary basis-size
      cpassert(SIZE(v_blocks, 1) == n .AND. SIZE(v_blocks, 2) == n)
      kterm = 0 ! init counter
 
      DO ipot = 1, npots
         pot_maxl = ipao_potentials(ipot)%maxl ! taken from central atom
 
         ! setup description of potential
         lb_min = 0
         lb_max = pot_maxl
         ncfgb = ncoset(lb_max) - ncoset(lb_min - 1)
         npgfb = 1 ! number of exponents
         nb = npgfb*ncfgb
         ALLOCATE (rpgfb(npgfb), zetb(npgfb))
 
         ! build block_V_full
         ALLOCATE (block_v_full(n, n, pot_maxl/2 + 1))
         block_v_full = 0.0_dp
 
         DO iset = 1, basis_set%nset
         DO jset = 1, iset
 
            ! setup iset
            la1_max = basis_set%lmax(iset)
            la1_min = basis_set%lmin(iset)
            npgfa1 = basis_set%npgf(iset)
            ncfga1 = ncoset(la1_max) - ncoset(la1_min - 1)
            na1 = npgfa1*ncfga1
            zeta1 => basis_set%zet(:, iset)
            rpgfa1 => basis_set%pgf_radius(:, iset)
 
            ! setup jset
            la2_max = basis_set%lmax(jset)
            la2_min = basis_set%lmin(jset)
            npgfa2 = basis_set%npgf(jset)
            ncfga2 = ncoset(la2_max) - ncoset(la2_min - 1)
            na2 = npgfa2*ncfga2
            zeta2 => basis_set%zet(:, jset)
            rpgfa2 => basis_set%pgf_radius(:, jset)
 
            ! radius of most diffuse basis-function
            rpgfa_max = max(maxval(rpgfa1), maxval(rpgfa2))
 
            ! allocate space for integrals
            ALLOCATE (saab(na1, na2, nb), saal(na1, na2, pot_maxl/2 + 1))
            saal = 0.0_dp
 
            ! loop over neighbors
            DO jatom = 1, natoms
               IF (jatom == iatom) cycle ! no self-interaction
               CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
               CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=jpao_potentials)
               IF (SIZE(jpao_potentials) /= npots) THEN
                  cpabort("Not all KINDs have the same number of PAO_POTENTIAL sections")
               END IF
 
               ! initialize exponents
               pot_weight = jpao_potentials(ipot)%weight ! taken from remote atom
               pot_beta = jpao_potentials(ipot)%beta ! taken from remote atom
               rpgfb(1) = jpao_potentials(ipot)%beta_radius ! taken from remote atom
               zetb(1) = pot_beta
 
               ! calculate direction
               ra = particle_set(iatom)%r
               rb = particle_set(jatom)%r
               rab = pbc(ra, rb, cell)
 
               ! distance screening
               tab = sqrt(sum(rab**2))
               IF (rpgfa_max + rpgfb(1) < tab) cycle
 
               ! calculate actual integrals
               saab = 0.0_dp
               CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
                                la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
                                lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
                                rab=rab, saab=saab)
 
               ! sum neighbor contributions according to remote atom's weight and normalization
               DO lpot = 0, pot_maxl, 2
                  norm2 = (2.0_dp*pot_beta)**(-0.5_dp - lpot)*gamma1(lpot)
                  ! sum potential terms: POW(x**2 + y**2 + z**2, lpot/2)
                  DO ic = ncoset(lpot - 1) + 1, ncoset(lpot)
                     coeff = multinomial(lpot/2, indco(:, ic)/2)
                     saal(:, :, lpot/2 + 1) = saal(:, :, lpot/2 + 1) + saab(:, :, ic)*coeff*pot_weight/sqrt(norm2)
                  END DO
               END DO
            END DO ! jatom
 
            ! find bounds of set-pair and setup transformation matrices
            sgfa1 = basis_set%first_sgf(1, iset)
            sgla1 = sgfa1 + basis_set%nsgf_set(iset) - 1
            sgfa2 = basis_set%first_sgf(1, jset)
            sgla2 = sgfa2 + basis_set%nsgf_set(jset) - 1
            t1 => basis_set%scon(1:na1, sgfa1:sgla1)
            t2 => basis_set%scon(1:na2, sgfa2:sgla2)
 
            ! transform into primary basis
            DO lpot = 0, pot_maxl, 2
               v12 => block_v_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1)
               v21 => block_v_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1)
               v12 = matmul(transpose(t1), matmul(saal(:, :, lpot/2 + 1), t2))
               v21 = transpose(v12)
            END DO
            DEALLOCATE (saab, saal)
         END DO ! jset
         END DO ! iset
         DEALLOCATE (rpgfb, zetb)
 
         ! block_V_full is ready -------------------------------------------------------------------
         ! split the full blocks into shell-pair sub-blocks
         DO iset = 1, basis_set%nset
         DO jset = 1, iset
         DO ishell = 1, basis_set%nshell(iset)
         DO jshell = 1, basis_set%nshell(jset)
            IF (basis_set%l(ishell, iset) == 0 .AND. basis_set%l(jshell, jset) == 0) cycle ! covered by central terms
            ishell_abs = sum(basis_set%nshell(1:iset - 1)) + ishell
            jshell_abs = sum(basis_set%nshell(1:jset - 1)) + jshell
            IF (ishell_abs < jshell_abs) cycle
 
            ! find bounds of shell-pair
            sgfa1 = basis_set%first_sgf(ishell, iset)
            sgla1 = basis_set%last_sgf(ishell, iset)
            sgfa2 = basis_set%first_sgf(jshell, jset)
            sgla2 = basis_set%last_sgf(jshell, jset)
 
            DO lpot = 0, pot_maxl, 2
               kterm = kterm + 1
               v_blocks(:, :, kterm) = 0.0_dp
               v_blocks(sgfa1:sgla1, sgfa2:sgla2, kterm) = block_v_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1)
               v_blocks(sgfa2:sgla2, sgfa1:sgla1, kterm) = block_v_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1)
            END DO ! lpot
         END DO ! jshell
         END DO ! ishell
         END DO ! jset
         END DO ! iset
         DEALLOCATE (block_v_full)
      END DO ! ipot
 
      ! terms on central atom ----------------------------------------------------------------------
 
      DO iset = 1, basis_set%nset
      DO jset = 1, iset
      DO ishell = 1, basis_set%nshell(iset)
      DO jshell = 1, basis_set%nshell(jset)
         IF (basis_set%l(ishell, iset) /= basis_set%l(jshell, jset)) cycle ! need quadratic block
         ishell_abs = sum(basis_set%nshell(1:iset - 1)) + ishell
         jshell_abs = sum(basis_set%nshell(1:jset - 1)) + jshell
         IF (ishell_abs < jshell_abs) cycle
         kterm = kterm + 1
         sgfa1 = basis_set%first_sgf(ishell, iset)
         sgla1 = basis_set%last_sgf(ishell, iset)
         sgfa2 = basis_set%first_sgf(jshell, jset)
         sgla2 = basis_set%last_sgf(jshell, jset)
         cpassert((sgla1 - sgfa1) == (sgla2 - sgfa2)) ! should be a quadratic block
         v_blocks(:, :, kterm) = 0.0_dp
         DO i = 1, sgla1 - sgfa1 + 1 ! set diagonal of sub-block
            v_blocks(sgfa1 - 1 + i, sgfa2 - 1 + i, kterm) = 1.0_dp
            v_blocks(sgfa2 - 1 + i, sgfa1 - 1 + i, kterm) = 1.0_dp
         END DO
         norm2 = sum(v_blocks(:, :, kterm)**2)
         v_blocks(:, :, kterm) = v_blocks(:, :, kterm)/sqrt(norm2) ! normalize
      END DO ! jshell
      END DO ! ishell
      END DO ! jset
      END DO ! iset
 
      cpassert(SIZE(v_blocks, 3) == kterm) ! ensure we generated all terms
 
      CALL timestop(handle)
   END SUBROUTINE linpot_rotinv_calc_terms
 
! **************************************************************************************************
!> \brief Calculate force contribution from rotinv parametrization
!> \param qs_env ...
!> \param iatom ...
!> \param M_blocks ...
!> \param forces ...
! **************************************************************************************************
   SUBROUTINE linpot_rotinv_calc_forces(qs_env, iatom, M_blocks, forces)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      INTEGER, INTENT(IN)                                :: iatom
      REAL(dp), DIMENSION(:, :, :), INTENT(IN)           :: m_blocks
      REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: forces
 
      CHARACTER(len=*), PARAMETER :: routinen = 'linpot_rotinv_calc_forces'
 
      INTEGER :: handle, i, ic, ikind, ipot, iset, ishell, ishell_abs, jatom, jkind, jset, jshell, &
         jshell_abs, kterm, la1_max, la1_min, la2_max, la2_min, lb_max, lb_min, lpot, n, na1, na2, &
         natoms, nb, ncfga1, ncfga2, ncfgb, npgfa1, npgfa2, npgfb, npots, nshells, pot_maxl, &
         sgfa1, sgfa2, sgla1, sgla2
      REAL(dp)                                           :: coeff, f, norm2, pot_beta, pot_weight, &
                                                            rpgfa_max, tab
      REAL(dp), DIMENSION(3)                             :: ra, rab, rb
      REAL(dp), DIMENSION(:), POINTER                    :: rpgfa1, rpgfa2, rpgfb, zeta1, zeta2, zetb
      REAL(dp), DIMENSION(:, :), POINTER                 :: block_d, t1, t2
      REAL(dp), DIMENSION(:, :, :), POINTER              :: block_m_full, dab
      REAL(dp), DIMENSION(:, :, :, :), POINTER           :: daab
      TYPE(cell_type), POINTER                           :: cell
      TYPE(gto_basis_set_type), POINTER                  :: basis_set
      TYPE(pao_potential_type), DIMENSION(:), POINTER    :: ipao_potentials, jpao_potentials
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, &
                      natom=natoms, &
                      cell=cell, &
                      particle_set=particle_set, &
                      qs_kind_set=qs_kind_set)
 
      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
      CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_potentials=ipao_potentials)
      npots = SIZE(ipao_potentials)
      nshells = sum(basis_set%nshell)
      n = basis_set%nsgf ! primary basis-size
      cpassert(SIZE(m_blocks, 1) == n .AND. SIZE(m_blocks, 2) == n)
      kterm = 0 ! init counter
      ALLOCATE (block_d(n, n))
 
      DO ipot = 1, npots
         pot_maxl = ipao_potentials(ipot)%maxl ! taken from central atom
 
         ! build block_M_full
         ALLOCATE (block_m_full(n, n, pot_maxl/2 + 1))
         block_m_full = 0.0_dp
         DO iset = 1, basis_set%nset
         DO jset = 1, iset
            DO ishell = 1, basis_set%nshell(iset)
            DO jshell = 1, basis_set%nshell(jset)
               IF (basis_set%l(ishell, iset) == 0 .AND. basis_set%l(jshell, jset) == 0) cycle ! covered by central terms
               ishell_abs = sum(basis_set%nshell(1:iset - 1)) + ishell
               jshell_abs = sum(basis_set%nshell(1:jset - 1)) + jshell
               IF (ishell_abs < jshell_abs) cycle
               ! find bounds of shell-pair
               sgfa1 = basis_set%first_sgf(ishell, iset)
               sgla1 = basis_set%last_sgf(ishell, iset)
               sgfa2 = basis_set%first_sgf(jshell, jset)
               sgla2 = basis_set%last_sgf(jshell, jset)
               DO lpot = 0, pot_maxl, 2
                  kterm = kterm + 1
                  block_m_full(sgfa1:sgla1, sgfa2:sgla2, lpot/2 + 1) = m_blocks(sgfa1:sgla1, sgfa2:sgla2, kterm)
                  block_m_full(sgfa2:sgla2, sgfa1:sgla1, lpot/2 + 1) = m_blocks(sgfa2:sgla2, sgfa1:sgla1, kterm)
               END DO ! lpot
            END DO ! jshell
            END DO ! ishell
         END DO ! jset
         END DO ! iset
 
         ! setup description of potential
         lb_min = 0
         lb_max = pot_maxl
         ncfgb = ncoset(lb_max) - ncoset(lb_min - 1)
         npgfb = 1 ! number of exponents
         nb = npgfb*ncfgb
         ALLOCATE (rpgfb(npgfb), zetb(npgfb))
 
         DO iset = 1, basis_set%nset
         DO jset = 1, iset
 
            ! setup iset
            la1_max = basis_set%lmax(iset)
            la1_min = basis_set%lmin(iset)
            npgfa1 = basis_set%npgf(iset)
            ncfga1 = ncoset(la1_max) - ncoset(la1_min - 1)
            na1 = npgfa1*ncfga1
            zeta1 => basis_set%zet(:, iset)
            rpgfa1 => basis_set%pgf_radius(:, iset)
 
            ! setup jset
            la2_max = basis_set%lmax(jset)
            la2_min = basis_set%lmin(jset)
            npgfa2 = basis_set%npgf(jset)
            ncfga2 = ncoset(la2_max) - ncoset(la2_min - 1)
            na2 = npgfa2*ncfga2
            zeta2 => basis_set%zet(:, jset)
            rpgfa2 => basis_set%pgf_radius(:, jset)
 
            ! radius of most diffuse basis-function
            rpgfa_max = max(maxval(rpgfa1), maxval(rpgfa2))
 
            ! find bounds of set-pair and setup transformation matrices
            sgfa1 = basis_set%first_sgf(1, iset)
            sgla1 = sgfa1 + basis_set%nsgf_set(iset) - 1
            sgfa2 = basis_set%first_sgf(1, jset)
            sgla2 = sgfa2 + basis_set%nsgf_set(jset) - 1
            t1 => basis_set%scon(1:na1, sgfa1:sgla1)
            t2 => basis_set%scon(1:na2, sgfa2:sgla2)
 
            ! allocate space for integrals
            ALLOCATE (daab(na1, na2, nb, 3), dab(na1, na2, 3))
 
            ! loop over neighbors
            DO jatom = 1, natoms
               IF (jatom == iatom) cycle ! no self-interaction
               CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
               CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=jpao_potentials)
               IF (SIZE(jpao_potentials) /= npots) THEN
                  cpabort("Not all KINDs have the same number of PAO_POTENTIAL sections")
               END IF
 
               ! initialize exponents
               pot_weight = jpao_potentials(ipot)%weight ! taken from remote atom
               pot_beta = jpao_potentials(ipot)%beta ! taken from remote atom
               rpgfb(1) = jpao_potentials(ipot)%beta_radius ! taken from remote atom
               zetb(1) = pot_beta
 
               ! calculate direction
               ra = particle_set(iatom)%r
               rb = particle_set(jatom)%r
               rab = pbc(ra, rb, cell)
 
               ! distance screening
               tab = sqrt(sum(rab**2))
               IF (rpgfa_max + rpgfb(1) < tab) cycle
 
               ! calculate actual integrals
               daab = 0.0_dp
               CALL overlap_aab(la1_max=la1_max, la1_min=la1_min, npgfa1=npgfa1, rpgfa1=rpgfa1, zeta1=zeta1, &
                                la2_max=la2_max, la2_min=la2_min, npgfa2=npgfa2, rpgfa2=rpgfa2, zeta2=zeta2, &
                                lb_max=lb_max, lb_min=lb_min, npgfb=npgfb, rpgfb=rpgfb, zetb=zetb, &
                                rab=rab, daab=daab)
 
               ! sum neighbor contributions according to remote atom's weight and normalization
               DO lpot = 0, pot_maxl, 2
                  ! sum potential terms: POW(x**2 + y**2 + z**2, lpot/2)
                  dab = 0.0_dp
                  DO ic = ncoset(lpot - 1) + 1, ncoset(lpot)
                     norm2 = (2.0_dp*pot_beta)**(-0.5_dp - lpot)*gamma1(lpot)
                     coeff = multinomial(lpot/2, indco(:, ic)/2)
                     dab = dab + coeff*daab(:, :, ic, :)*pot_weight/sqrt(norm2)
                  END DO
                  DO i = 1, 3
                     ! transform into primary basis
                     block_d = 0.0_dp
                     block_d(sgfa1:sgla1, sgfa2:sgla2) = matmul(transpose(t1), matmul(dab(:, :, i), t2))
                     block_d(sgfa2:sgla2, sgfa1:sgla1) = transpose(block_d(sgfa1:sgla1, sgfa2:sgla2))
                     ! calculate and add forces
                     f = sum(block_m_full(:, :, lpot/2 + 1)*block_d)
                     forces(iatom, i) = forces(iatom, i) - f
                     forces(jatom, i) = forces(jatom, i) + f
                  END DO
               END DO ! lpot
            END DO ! jatom
            DEALLOCATE (dab, daab)
         END DO ! jset
         END DO ! iset
         DEALLOCATE (rpgfb, zetb, block_m_full)
      END DO ! ipot
      DEALLOCATE (block_d)
 
      CALL timestop(handle)
   END SUBROUTINE linpot_rotinv_calc_forces
 
END MODULE pao_linpot_rotinv