Among several possible approaches we can view this problem as a
straightforward application of the Polya Enumeration Theorem and very
similar to the problem discussed at this MSE
link.
Recall the recurrence by Lovasz for the cycle index $Z(S_n)$ of
the multiset operator $\mathfrak{M}_{=n}$ on $n$ slots, which is
$$Z(S_n) = \frac{1}{n} \sum_{l=1}^n a_l Z(S_{n-l})
\quad\text{where}\quad
Z(S_0) = 1.$$
Suppose the prime factorization of $n$ is given by
$$n=\prod_p p^v.$$
Applying PET it now follows almost by inspection that the desired
count of multiplicative partitions into $k$ factors is given by the
term
$$F(n, k) = \left[\prod_p X_p^v\right]
Z(S_k)\left(-1 + \prod_p \frac{1}{1-X_p}\right)$$
where the square bracket denotes coefficient extraction
of formal power series.
We are interested in
$$G(n) = \sum_{k=1}^{\Omega(n)} F(n, k).$$
Now recall the generating function of the cycle index of the symmetric
group which is
$$H(z) = \sum_{q\ge 0} Z(S_q) z^q =
\exp\left(a_1z + a_2\frac{z^2}{2} + a_3\frac{z^3}{3}+\cdots\right)
= \exp\left(\sum_{l\ge 1} a_l \frac{z^l}{l}\right).$$
This gives for $G(n)$
$$G(n) = \left[\prod_p X_p^v\right]
\exp\left(\sum_{l\ge 1} \frac{1}{l}
\left(-1 + \prod_p \frac{1}{1-X_p^l}\right)\right).$$
Recalling $v_p(n)$ which is the exponent of the maximum power of $p$
that divides $n$ we finally obtain
$$G(n) = \left[\prod_p X_p^v\right]
\exp\left(\sum_{l=1}^{\max_p v_p(n)} \frac{1}{l}
\left(-1 + \prod_p \frac{1}{1-X_p^l}\right)\right).$$
This formula has a certain aesthetic value. Implemented in Maple it
will produce the following sequence:
$$1, 1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 4, 1, 2, 2, 5, 1, 4, 1, 4,
\\ 2, 2, 1, 7, 2, 2, 3, 4, 1, 5, 1, 7, 2, 2, 2, 9, 1, 2, 2, 7,
\\ 1, 5, 1, 4, 4, 2, 1, 12, \ldots$$
which points us to OEIS A001055 where
additional information awaits, including simple recurrences for
practical computation of this function.
The Maple code for the above formula is as follows (here we have
replaced the prime number index on the variable $X$ by a positional
index on the variable $Y$):
with(numtheory);
facts :=
proc(n)
option remember;
local f, gf, p, res;
f := op(2, ifactors(n));
gf := exp(add(1/l*
(-1+mul(1/(1-Y[p]^l), p=1..nops(f))),
l=1..max(map(p->p[2], f))));
res := gf;
for p to nops(f) do
res := coeftayl(res, Y[p]=0, f[p][2]);
od;
res;
end;