def _fcq_to_fcr(
        self, fcq: dict[tuple[float, float, float], np.ndarray]
    ) -> np.ndarray:
        """
        Fourier transform the force constants on the full q-grid to real
        space.
        """
        all_q = np.array(list(fcq.keys()))
        print(f"Found {len(all_q)} q-points in the full BZ")
        dynmat = np.array([self._fcq_to_dynmat(fcq[tuple(q)], q) for q in all_q])
        d2f = DynmatToForceConstants(
            self.sc_gen.primitive,
            self.sc_gen.supercell,
            is_full_fc=True,
        )
        d2f.commensurate_points = all_q
        d2f.dynamical_matrices = dynmat
        print("Running DynmatToForceConstants...")
        d2f.run()
        return d2f.force_constants

  nq = 0
  do isym = 1, nsym
    ism1 = invs(isym)
    do i = 1, 3
      raq(i) = s(i, 1, ism1)*aq(1) &
               + s(i, 2, ism1)*aq(2) &
               + s(i, 3, ism1)*aq(3)
    end do
    do i = 1, 3
      sxq(i, 48) = bg(i, 1)*raq(1) &
                   + bg(i, 2)*raq(2) &
                   + bg(i, 3)*raq(3)
    end do
    do iq = 1, nq
      if (eqvect(raq, saq(1, iq), zero)) then
        isq(isym) = iq
        nsq(iq) = nsq(iq) + 1
      end if
    end do
    if (isq(isym) == 0) then
      nq = nq + 1
      nsq(nq) = 1
      isq(isym) = nq
      saq(:, nq) = raq(:)
      do i = 1, 3
        sxq(i, nq) = bg(i, 1)*saq(1, nq) &
                     + bg(i, 2)*saq(2, nq) &
                     + bg(i, 3)*saq(3, nq)
      end do
    end if
  end do

function ℒ(r, p)
    rᵢⱼ = pairwise(Euclidean(), r')
    mask = .!I(p.N_atoms) # Mask to avoid r_{i=j} term

    # Subtraction accounts for the r_{i=j} term
    X = sum(k.(rᵢⱼ, p.rμ, p.rσ), dims=1) .- k(0, p.rμ, p.rσ) - p.m'
    ℒC = sum(abs2, X) # Coordination number
    ℒmin = sum(1 .- 1 ./ (1 .+ exp.(-p.γ .* (rᵢⱼ[mask] .- p.rmin) ./ p.a))) # Minimum distance
    #ℒmax = sum(exp.((rᵢⱼ[mask] .- p.rmax) ./ p.a) ./ (1 .+ exp.(-p.β .* (rᵢⱼ[mask] .- p.rmax) / p.a))) # Maximum distance
    ℒmax = sum(exp.((rᵢⱼ[mask] .- p.rmax)) ./ (1 .+ exp.(-p.β .* (rᵢⱼ[mask]- p.rmax)))) # Maximum distance
    @show ℒC, ℒmin, ℒmax
    return p.λ₁ * ℒC + p.λ₂ * ℒmin + p.λ₃ * ℒmax
end