As a matter of fact, it isn't hard to construct a multiplicative sequence $a_n$ such that $f(z)$ is an entire function without zeroes. Unfortunately, it is completely useless for the questions that you brought up as "motivation".
Here is the construction.
Claim 1: Let $\lambda_j\in [0,1]$ ($j=0,\dots,M$). Assume that $|a_j|\le 1$ and $\sum_ja_j\lambda_j^p=0$ for all $0\le p\le P$. Then, if $P>2eT$, we have $\left|\sum_j a_je^{\lambda_j z}\right|\le (M+1)(eT/P)^P$
Proof: Taylor decomposition and a straightforward tail estimate.
Claim 2: Let $P$ be large enough. Let $\Delta>0$, $M>P^3$ and $(M+1)\Delta<1$. Let $I_j$ ($0\le j\le M$) be $M+1$ adjacent intervals of length about $\Delta$ each arranged in the increasing order such that $I_0$ contains $0=\lambda_0$. Suppose that we choose one $\lambda_j$ in each interval $I_j$ with $j\ge 1$. Then for every $|a_0|\le 1$, there exist $a_j\in\mathbb C$ such that $|a_j|\le 1$ and $\sum_{j\ge 0} a_j\lambda_j^p=0$ for $0\le p\le P$.
Proof: By duality, we can restate it as the claim that $\sum_{j\ge 1}|Q(\lambda_j)|\ge |Q(0)|$ for every polynomial $Q$ of degree $P$. Now, let $I$ be the union of $I_j$. It is an interval of length about $M\Delta$, so, by Markov's inequality, $|Q'|\le CP^2(M\Delta)^{-1} A\le CP^{-1}\Delta^{-1}A$ where $A=\max_I |Q|\ge |Q(0)|$. But then on the 5 intervals $I_j$ closest to the point where the maximum is attained, we have $|Q|\ge A-5\Delta CP^{-1}\Delta^{-1}A\ge A/2$. The rest is obvious.
Claim 3: Suppose that $a_0$ is fixed and $a_j$ ($j\ge 1$) satisfy $\sum_{j\ge 0} a_j\lambda_j^p=0$ for $0\le p\le P$. Then we can change $a_j$ with $j\ge 1$ so that all but $P+1$ of them are exactly $1$ in absolute value and the identities still hold.
Proof: As long as we have more than $P+1$ small $a_j$, we have a line of solutions of our set of $P+1$ linear equations. Moving along this line we can make one of small $a_j$ large. Repeating this as long as possible, we get the claim.
Now it is time to recall that the logarithm of the function $f(z)$ is given by
$$
L(z)=\sum_{n\in\Lambda}a_ne^{-z\log n}
$$
where $\Lambda$ is the set of primes and prime powers and $a_n=m^{-1}a_p^m$ if $n=p^m$. We are free to choose $a_p$ for prime $p$ in any way we want but the rest $a_n$ will be uniquely determined then. The key point is that we have much more primes than prime powers for unit length.
So, split big positive numbers into intervals from $u$ to $e^\Delta u$ where $\Delta$ is a slowly decaying function of $u$ (we'll specify it later). Formally we define the sequence $u_k$ by $u_0=$something large, $u_{k+1}=e^{\Delta(u_k)}u_k$ but to put all those backward apostrophes around formulae is too big headache, so I'll drop all indices. Choose also some slowly growing functions $M=M(u)$ and $P=P(u)$ to be specified later as well.
We need a few things:
1) Each interval should contain many primes. Since the classical prime number theorem has the error term $ue^{-c\sqrt{\log u}}$, this calls for $\Delta=\exp\{-\log^{\frac 13} u\}$
Then we still have at least $u^{4/5}$ primes in each interval (all we need is to beat $u^{1/2}$ with some margin).
2) We should have $M\Delta\ll 1$, $M>P^3$, and $u\left(\frac{eT}P\right)^P\le (2u)^{-T-3}$ for any fixed $T>0$ and all sufficiently large $u$. This can be easily achieved by choosing $P=\log^2 u$ and $M=\log^6 u$.
3) At last, we'll need $M(P+\sqrt u)\ll u^{4/5}$, which is true for our choice.
Now it is time to run the main inductive construction. Suppose that $a_n$ are already chosen for all $n$ in the intervals up to $(u,e^{\Delta}u)$ and we still have almost all primes in the intervals following the current interval free (we'll check this condition in the end). We want to assign $a_p$ for all $p$ in our interval for which the choice hasn't been made yet or was made badly. We start with looking at all $a_p$ that are not assigned yet or assigned in a lame way, i.e., less than one in absolute value. Claim 3 (actually a small modification of it) allows us to upgrade all of them but $P+1$ to good ones (having absolute value $1$) at the expense of adding an entire function that in the disk of radius $T$ is bounded by $(2u)^{T} u\left(\frac{eT}P\right)^P\le u^{-3}$ to $L(z)$. Now we are left with at most $\sqrt u$ powers of primes and $P+1$ lame primes to take care of. We need the prime powers participate in small sums as they are and we need the small coefficients to be complemented by something participating in small sums too. For each of them, we choose $M$ still free primes in the next $M$ intervals (one in each) and apply Claim 2 to make a (lame) assignment so that the corresponding sum is again bounded by $u^{-3}$ in the disk of radius $T$. We have at most $u$ such sums, so the total addition will be at most $u^{-2}$. This will finish the interval off. Now it remains to notice that we used only about $\sqrt u+P$ free primes in each next interval and went only $M$ intervals ahead. This means that in each interval only $M(\sqrt u+P)$ free primes will ever be used for compensating the previous intervals, so we'll never run out of free primes. Also, the sum of the blocks we constructed will converge to an entire function. At last, when $\Re z>1$, we can change the order of summation and exponentiate finally getting the Dirichlet series representation that we need.
The end.
