new_lowrank_c6 Subroutine

subroutine new_lowrank_c6(self, ref, c6ref, config)

   !> Instance of the separable representation
   type(d3_lowrank_c6), intent(out) :: self

   !> Number of reference systems for each species
   integer, intent(in) :: ref(:)

   !> Reference C6 coefficients
   real(wp), intent(in) :: c6ref(:, :, :, :)

   !> Setup of the separable representation
   type(d3_lowrank_config), intent(in) :: config

   integer :: nid, mref, ndim, isp, iref, ii, jj, il, nrank
   integer, allocatable :: spmap(:), refmap(:)
   real(wp), allocatable :: amat(:, :), eval(:), evec(:, :), res(:, :)

   nid = size(ref)
   mref = maxval(ref)
   ndim = sum(ref)

   allocate(spmap(ndim), refmap(ndim))
   ii = 0
   do isp = 1, nid
      do iref = 1, ref(isp)
         ii = ii + 1
         spmap(ii) = isp
         refmap(ii) = iref
      end do
   end do

   allocate(amat(ndim, ndim))
   do ii = 1, ndim
      do jj = 1, ndim
         amat(jj, ii) = c6ref(refmap(jj), refmap(ii), spmap(jj), spmap(ii))
      end do
   end do

   allocate(eval(ndim), evec(ndim, ndim))
   call symmetric_eigendecomposition(amat, eval, evec)

   allocate(res(ndim, ndim), source=amat)
   nrank = ndim
   do il = 1, ndim
      do ii = 1, ndim
         res(:, ii) = res(:, ii) - eval(il) * evec(:, il) * evec(ii, il)
      end do
      self%error = max_relative_error(res, amat)
      if (config%rank > 0) then
         if (il >= min(config%rank, ndim)) then
            nrank = il
            exit
         end if
      else if (self%error <= config%tolerance) then
         nrank = il
         exit
      end if
   end do

   self%rank = nrank
   self%kcut = config%kcut
   self%mesh = config%mesh
   self%lambda = eval(:nrank)
   allocate(self%vec(mref, nid, nrank), source=0.0_wp)
   do il = 1, nrank
      do ii = 1, ndim
         self%vec(refmap(ii), spmap(ii), il) = evec(ii, il)
      end do
   end do

end subroutine new_lowrank_c6