This is a straightforward application of the Polya Enumeration
Theorem. We treat the problem of subsets with $n$ elements of the set
$\{1,2,\ldots, q\}$ whose sum is divisible by $k.$ Suppose $Z(P_n)$ is
the cycle index of the set operator $\mathfrak{P}_{=n}$ given by the
recurrence by Lovasz which is
$$Z(P_n) = \frac{1}{n} \sum_{l=1}^n (-1)^{l+1} a_l Z(P_{n-l})
\quad\text{where}\quad
Z(P_0) = 1.$$
We obtain by PET the following formula for the OGF of ordinary sets
$$Z(P_n)\left(w+w^2+\cdots+w^q\right)
= Z(P_n)\left(\sum_{m=1}^q w^m\right).$$
With $\rho$ a root of unity namely
$$\rho = \exp(2\pi i/k)$$
we get for the desired count the value
$$\frac{1}{k}
\left.\sum_{p=0}^{k-1}
Z(P_n)\left(\sum_{m=1}^q w^m\right)\right|_{w=\rho^p}.$$
We will compute the value for $p=0$ separately and to do this recall
the OGF of the set operator $\mathfrak{P}_{=n}$ which is
$$Z(P_n) = [z^n]
\exp\left(a_1 z - a_2 \frac{z^2}{2}
+ a_3 \frac{z^3}{3}
- a_4 \frac{z^4}{4}
+\cdots \right).$$
or
$$Z(P_n) =
[z^n] \exp\left(\sum_{r\ge 1} (-1)^{r+1} a_r \frac{z^r}{r}\right).$$
On substituting this into our formula we get
$$\frac{1}{k}
\sum_{p=0}^{k-1}
\left. [z^n]
\exp\left(\sum_{r\ge 1} (-1)^{r+1} \frac{z^r}{r}
\sum_{m=1}^q w^{rm} \right) \right|_{w=\rho^p}.$$
When $p=0$ we obtain
$$\frac{1}{k}
[z^n] \exp\left(q\sum_{r\ge 1} (-1)^{r+1} \frac{z^r}{r}\right)
= \frac{1}{k}
[z^n] \exp\left(q \log (1+z)\right)
\\ = \frac{1}{k} [z^n] (1+z)^q
= \frac{1}{k} {q\choose n}.$$
We switch to algorithmics for the remainder of this discussion.
In treating the case $p\ge 1$ we make the following observation.
When substituting into the terms of the cycle index those $a_r$ from
the product where $pr$ is a multiple of $k$ produce the value $q$
while the remaining $a_r$ create a sequence of period $k$ that depends
only on the remainder $b$ when $q$ is divided by $k$ where we take
$1\le b\le k.$
This yields an algorithm where we iterate over the cycle index,
extract eventual powers of $q$ from the terms and interpolate the rest
in terms of $b.$ The algorithm can be used to compute formulae for
fixed combinations of $n$ and $k$ like the ones at this MSE
link automatically.
We obtain for $(n,k) = (3,3)$
$$1/18\,{q}^{3}+1/3\,{b}^{2}-1/6\,{q}^{2}-{\frac {11\,b}{9}}+q/3+2/3$$
and sure enough comparing it to the link these are the right values.
Supposing now that we are interested in divisibility by five of
three-element subsets i.e. the pair $(n,k) = (3,5)$ we find
$$1/12\,{b}^{4}-{\frac {13\,{b}^{3}}{15}}+1/30\,{q}^{3}+{\frac {
181\,{b}^{2}}{60}}-1/10\,{q}^{2}-{\frac {127\,b}{30}}+q/15+2$$
which gives the sequence (starting at $q=3$)
$$0, 0, 2, 4, 7, 11, 16, 24, 33, 44, 57, 72, 91, \ldots$$
For the pair $(11,5)$ we obtain
$${\frac {{q}^{11}}{199584000}}-{\frac {{q}^{10}}{3628800}}+{
\frac {{q}^{9}}{151200}}-{\frac {11\,{q}^{8}}{120960}}+{\frac {
683\,{q}^{7}}{864000}}+{\frac {{b}^{4}{q}^{2}}{1200}}\\-{\frac {
781\,{q}^{6}}{172800}}-{\frac {{b}^{4}q}{80}}-{\frac {{b}^{3}{q
}^{2}}{120}}+{\frac {31063\,{q}^{5}}{1814400}}+1/24\,{b}^{4}+1/
8\,{b}^{3}q+{\frac {7\,{b}^{2}{q}^{2}}{240}}\\-{\frac {1529\,{q}^
{4}}{36288}}-{\frac {631\,{b}^{3}}{1500}}-{\frac {859\,{b}^{2}q
}{2000}}-{\frac {137\,b{q}^{2}}{3000}}+{\frac {16103\,{q}^{3}}{
252000}}+{\frac {863\,{b}^{2}}{600}}\\+{\frac {129\,bq}{200}}-{
\frac {419\,{q}^{2}}{12600}}-{\frac {25\,b}{12}}-{\frac {31\,q
}{110}}+1$$
which gives the sequence (starting at $q=11$)
$$0, 2, 15, 72, 273, 873, 2474, 6363, 15114, 33592, \ldots.$$
Another interesting pair is $(3,6)$ which gives
$$-{\frac {{b}^{5}q}{90}}-{\frac {{b}^{5}}{90}}+{\frac {7\,{b}^{4
}q}{36}}+1/4\,{b}^{4}-{\frac {23\,{b}^{3}q}{18}}-2\,{b}^{3}+{
\frac {35\,{b}^{2}q}{9}}\\+1/36\,{q}^{3}+7\,{b}^{2}-{\frac {242\,
bq}{45}}-1/12\,{q}^{2}-{\frac {313\,b}{30}}+{\frac {17\,q}{6}}+5$$
and $(4,7)$ which produces
$${\frac {{b}^{6}}{360}}-{\frac {3\,{b}^{5}}{40}}+{\frac {389\,{b
}^{4}}{504}}+{\frac {{q}^{4}}{168}}-{\frac {215\,{b}^{3}}{56}}-
1/28\,{q}^{3}\\+{\frac {3041\,{b}^{2}}{315}}+{\frac {11\,{q}^{2}
}{168}}-{\frac {403\,b}{35}}-q/28+5.$$
The Maple code for this including a total enumeration routine for
verification and some code to prettify the formulae for $k$ small was
as follows.
with(combinat);
pet_cycleind_set :=
proc(n)
option remember;
if n=0 then return 1; fi;
expand(1/n*add((-1)^(l+1)*a[l]*pet_cycleind_set(n-l), l=1..n));
end;
pet_flatten_term :=
proc(varp)
local terml, d, cf, v;
terml := [];
cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), seq(v, k=1..d)];
cf := cf/v^d;
od;
[cf, terml];
end;
V :=
proc(q, n, k)
local comb, res;
option remember;
res := 0;
comb := firstcomb(q, n);
while type(comb, set) do
if add(p, p in comb) mod k = 0 then
res := res + 1;
fi;
comb := nextcomb(comb, q);
od;
res;
end;
CF :=
proc(n, k)
local term, rho, res, p, flat, ixvar, rmd, m,
rep, qpow, w, ex, vals, val;
option remember;
res := 0; rho := exp(2*Pi*I/k);
for term in pet_cycleind_set(n) do
flat := pet_flatten_term(term);
for p to k-1 do
vals := []; rep := []; qpow := 0;
for ixvar in flat[2] do
ex := p*op(1, ixvar);
if ex mod k = 0 then
qpow := qpow + 1;
else
rep := [op(rep), ixvar];
fi;
od;
for rmd to k do
val := 1;
for ixvar in rep do
ex := p*op(1, ixvar);
val := val * add(w^(ex*m), m=1..rmd);
od;
vals :=
[op(vals), subs(w=rho, expand(val))];
od;
res := res + flat[1]*q^qpow *
interp([seq(b, b=1..k)], vals, b);
od;
od;
binomial(q,n)/k + res/k;
end;
CFsimp :=
proc(n, k)
local form, res, term, lcf, cnst;
option remember;
form := collect(expand(CF(n,k)), {b,q}, `distributed`);
cnst := coeff(coeff(form, b, 0), q, 0);
res := 0;
for term in form-cnst do
lcf := lcoeff(term);
res := res +
simplify(Re(lcf))*term/lcf;
od;
res + simplify(cnst);
end;
F :=
proc(qv, n, k)
local bv;
if qv mod k = 0 then
bv := k;
else
bv := qv mod k;
fi;
subs({q=qv, b=bv}, CFsimp(n,k));
end;
The reader is invited to contribute a better simplification routine
making more effective use of the mathematical givens. The Maple code
should be considered betaware.
Remark Sat Jan 23 2016. I present one of the special cases
where radical simplification is possible. Start from the formula
$$\frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
\left. [z^n]
\exp\left(\sum_{r\ge 1} (-1)^{r+1} \frac{z^r}{r}
\sum_{m=1}^q w^{rm} \right) \right|_{w=\rho^p}.$$
Now suppose that $q$ is a multiple of $k$ and $k$ is an odd
prime. Observe that the sum $$\sum_{m=1}^q w^{rm}$$ is equal to
$q/k\times k = q$ if $pr$ is a multiple of $k$ and zero otherwise. But
$pr$ can only be a multiple of $k$ if $r$ is a multiple of $k.$ This
yields
$$\frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
\left. [z^n]
\exp\left(\sum_{r\ge 1} (-1)^{kr+1} \frac{z^{kr}}{kr}
\sum_{m=1}^q w^{krm} \right) \right|_{w=\rho^p}
\\ = \frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
\left. [z^n]
\exp\left(\sum_{r\ge 1} (-1)^{kr+1} \frac{z^{kr}}{kr}
\frac{q}{k}\times k \right) \right|_{w=\rho^p}
\\ = \frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
[z^n]
\exp\left(\frac{q}{k}
\sum_{r\ge 1} (-1)^{kr+1} \frac{z^{kr}}{r}\right)
\\ = \frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
[z^n]
\exp\left(-\frac{q}{k}
\sum_{r\ge 1} \frac{(-z)^{kr}}{r}\right)
\\ = \frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
[z^n]
\exp\left(-\frac{q}{k} \log\frac{1}{1-(-z)^k}\right)
\\ = \frac{1}{k} {q\choose n} +
\frac{1}{k}
\sum_{p=1}^{k-1}
[z^n] (1+z^k)^{q/k}.$$
Therefore if $n$ is coprime to $k$ we obtain
$$\frac{1}{k} {q\choose n}$$
and if it is a multiple of $k$
$$\frac{1}{k} {q\choose n} +
\frac{k-1}{k} {q/k\choose n/k}.$$