grd2_tdhf_compute_data_t_get_density Subroutine

  subroutine grd2_tdhf_compute_data_t_get_density(this, basis, id, dab, dabmax)

    implicit none

    class(grd2_tdhf_compute_data_t), target, intent(inout) :: this
    type(basis_set), intent(in) :: basis
    integer, intent(in) :: id(4)
    real(kind=dp), target, intent(out) :: dab(*)
    real(kind=dp), intent(out) :: dabmax

    real(kind=dp) :: coul, coul2, coul3, exch, exch2, exch3, exch4
    real(kind=dp) :: df1
    real(kind=dp) :: coulfact, xcfact
    integer :: i_, j_, k_, l_
    integer :: i, j, k, l
    integer :: loc(4)
    integer :: nbf(4)
    real(kind=dp), pointer :: ab(:,:,:,:)

    coulfact = 4*this%coulscale
    xcfact = this%hfscale
    dabmax = 0
    loc = basis%ao_offset(id)-1

    nbf = basis%naos(id)

    ab(1:nbf(4),1:nbf(3),1:nbf(2),1:nbf(1)) => dab(1:product(nbf))

    associate( p => this%p2(:,:,1) &
             , d => this%d2(:,:,1) &
             , xpy => this%xpy2(:,:,1) &
             , xmy => this%xmy2(:,:,1) &
             )

    do i_ = 1, nbf(1)
      i = loc(1) + i_

      do j_ = 1, nbf(2)
        j = loc(2) + j_

        do k_ = 1, nbf(3)
          k = loc(3) + k_

          do l_ = 1, nbf(4)
            l = loc(4) + l_

            coul = &
            +   p(j,i) * d(l,k) &
            +   p(l,k) * d(j,i)

            coul2 = &
            + xpy(i,j) * xpy(k,l)

            coul3 = d(i,j) *   d(k,l)

            coul = 2*coul + 8*coul2 + coul3

            df1 = coulfact*coul

            if (xcfact/=0.0_dp) then
              exch = &
              +   p(i,k) *   d(l,j) &
              +   p(j,k) *   d(l,i) &
              +   p(l,i) *   d(k,j) &
              +   p(l,j) *   d(k,i)

              exch2 = &
              + xpy(k,i) * xpy(l,j) &
              + xpy(l,i) * xpy(k,j)

              exch3 = &
              + (xmy(i,l)-xmy(l,i)) * (xmy(j,k)-xmy(k,j)) &
              + (xmy(j,l)-xmy(l,j)) * (xmy(i,k)-xmy(k,i))

              exch4 =    d(i,k) *   d(j,l) &
                     +   d(i,l) *   d(j,k)

              exch = 2*exch + 8*exch2 + 2*exch3 + exch4

              df1 = df1 - xcfact*exch
            end if
            dabmax = max(dabmax, abs(df1))
            ab(l_,k_,j_,i_) = df1*product(basis%bfnrm([i,j,k,l]))
          end do
        end do
      end do
    end do

    end associate

  end subroutine grd2_tdhf_compute_data_t_get_density