Maximum value $c$ s.t. $\exists$ a subset $S$ of $\{z_1,z_2,\ldots,z_n\}$ s.t. $\left|\sum_{z\in S}z\right|\geq c$ ($\sum_{i=1}^{n}|z_i|=1$).
A bit of intuition to start: Our basic idea here will to be to think of our complex numbers as vectors in $\mathbb{R}^2$ and work with projections instead of lengths, mainly because projections work so nicely with respect to sums: The projection of the sum is the sum of the projections. So instead of trying to create a sum with large norm, we'll pick a direction and try to create a sum with large projection/component along that vector.
Suppose that I gave you a direction $v$ in advance, and all you wanted to do was maximize the component of your sum along $v$. Then choosing your subset would be easy: You include the vectors that point "towards" $v$ (have positive inner product), and exclude the ones that point "away" from $v$ (have negative inner product).
But since I'm the one picking $v$, you may be completely out of luck -- maybe I'll give you a $v$ that points away from all of your vectors, and the best component you could get is $0$. We want to pick the "right" $v$ in some sense -- one that gives us a large projection. But it's hard to do that without knowing the vectors in advance. And furthermore, Han de Bruijn's answer suggests that the extreme case comes when everything's symmetrical and we could just pick any direction whatsoever at random.
What this suggests is that we should pick a random direction, or to put it differently, show that the average $v$ gives us a large projection. Now for the actual calculation ...
Let $v$ be an arbitrary unit vector. As suggested in the intuition above, let $S_v$ denote those indices $i$ for which $z_i$ has a positive inner product with $v$, and define $$x_v := \sum_{i \in S_v} z_i.$$ We can bound $|x_v|$ from below by its component along $v$: $$|x_v| \geq \langle x_v, v\rangle = \sum_{i \in S_v} \langle z_i, v \rangle = \sum_{i \in S_v} |z_i| \cos(\theta_{i,v}),$$ where $\theta_{i,v}$ is the angle between $z_i$ and $v$. We can rewrite this as $$|x_v| \geq \sum_{i=1}^n |z_i| \max\{0, \cos(\theta_{i,v})\},$$ since by definition the extra terms we're including are all just $0$.
Now suppose I were to pick a $v$ uniformly at random from the unit circle. Then the angle $\theta_{i,v}$ would be uniform on $[0,2 \pi)$, and the expected value of the $i^{th}$ term on the right hand side would be $$\frac{1}{2 \pi} \int_0^{2\pi} |z_i| \max\{0, \cos \theta\} \, d\theta = \frac{1}{\pi} |z_i|.$$ Adding up over all $i$, we have $$E(|x_v|) \geq \frac{1}{\pi} \sum_{i=1}^n |z_i| = \frac{1}{\pi}$$ Here $E(|x_v|)$ denotes the expected value (or average) of $|x_v|$ taking over a randomly chosen $v$.
The key point here (what's sometimes termed "Erdős magic" after Paul Erdős): If we have a collection of vectors where the average length is $\frac{1}{\pi}$, that means there must be a vector in that collection whose length is at least $\frac{1}{\pi}$, so we win. That $\frac{1}{\pi}$ is the best possible can be shown using the point set from Han de Bruijn's answer.
This problem is a sort of old chestnut, and the argument here is not my own. But I don't actually know the original source for it. I'd love it if anyone who knows a bit more of its history could comment.
Special case first. It is supposed that the maximum value $\,c\,$ is achieved for evenly distributed $\{z_0,z_1,\ldots,z_{n-1}\}$ .
This means that
$$
z_k = \frac{1}{n} e^{i\cdot k\,2\pi/n} \quad \Longrightarrow \quad |z_0|+|z_1|+\cdots +|z_{n-1}|=1
$$
Easy results are obtained (by hand) for $n=1,2,3,4$:
$$
\begin{cases}
n=1 \quad \Longrightarrow \quad S = \{1\} & \mbox{and} & \left|\sum_{z\in S} z\right|\geq 1 > 1/\pi \\
n=2 \quad \Longrightarrow \quad S = \{1/2\} & \mbox{and} & \left|\sum_{z\in S} z\right|\geq 1/2 > 1/\pi \\
n=3 \quad \Longrightarrow \quad S = \{1/3,(-1-i\sqrt{3})/6\} & \mbox{and} & \left|\sum_{z\in S} z\right|\geq 1/3 > 1/\pi\\
n=4 \quad \Longrightarrow \quad S = \{1/4,i/4\} & \mbox{and} & \left|\sum_{z\in S} z\right|\geq \sqrt{2}/4 > 1/\pi
\end{cases}
$$
This already improves the bounds $1/6$ and $1/4$ as given in the question (at best for $n=3$).
Higher order results are obtained with help of a computer program. The lines giving $n$ , $\left|\sum_{z\in S} z\right|$
and $1/\pi$ are alternating with lines giving the indices $\,k\,$ of the terms $\,z_k\,$ in the sum $\,\left|\sum_{z\in S} z\right|$ .
It is shown that the sums indeed seem to converge to the conjectured value of $1/\pi$ . And there is a pattern in the subsets that do the job.
1 1.00000000000000E+0000 > 3.18309886183791E-0001
0
2 5.00000000000000E-0001 > 3.18309886183791E-0001
0
3 3.33333333333333E-0001 > 3.18309886183791E-0001
0 2
4 3.53553390593274E-0001 > 3.18309886183791E-0001
0 1
5 3.23606797749979E-0001 > 3.18309886183791E-0001
0 1 2
6 3.33333333333333E-0001 > 3.18309886183791E-0001
0 4 5
7 3.20997086245352E-0001 > 3.18309886183791E-0001
0 1 2
8 3.26640741219094E-0001 > 3.18309886183791E-0001
0 1 2 3
9 3.19931693507980E-0001 > 3.18309886183791E-0001
5 6 7 8
10 3.23606797749979E-0001 > 3.18309886183791E-0001
3 4 5 6 7
11 3.19394281060558E-0001 > 3.18309886183791E-0001
4 5 6 7 8 9
12 3.21975275429689E-0001 > 3.18309886183791E-0001
1 2 3 4 5 6
13 3.19085761944568E-0001 > 3.18309886183791E-0001
3 4 5 6 7 8
14 3.20997086245352E-0001 > 3.18309886183791E-0001
0 1 2 3 4 5 6
15 3.18892407783521E-0001 > 3.18309886183791E-0001
0 1 2 3 11 12 13 14
16 3.20364430967688E-0001 > 3.18309886183791E-0001
0 1 2 3 4 13 14 15
17 3.18763277866454E-0001 > 3.18309886183791E-0001
5 6 7 8 9 10 11 12
18 3.19931693507980E-0001 > 3.18309886183791E-0001
2 3 4 5 6 7 8 9 10
19 3.18672778564237E-0001 > 3.18309886183791E-0001
2 3 4 5 6 7 8 9 10
20 3.19622661074983E-0001 > 3.18309886183791E-0001
3 4 5 6 7 8 9 10 11 12
21 3.18606904753685E-0001 > 3.18309886183791E-0001
0 1 2 3 15 16 17 18 19 20
22 3.19394281060558E-0001 > 3.18309886183791E-0001
0 1 2 3 4 5 6 7 8 9 10
23 3.18557468338846E-0001 > 3.18309886183791E-0001
7 8 9 10 11 12 13 14 15 16 17
24 3.19220732314183E-0001 > 3.18309886183791E-0001
4 5 6 7 8 9 10 11 12 13 14 15
Two snapshots should clarify the pattern in the subsets that do the job:
So it seems that, without too much loss of generality, we can conjecture that:
$$
\left|\sum_{z\in S} z\right| = \left|\sum_{k=0}^{n/2-1} \frac{1}{n} e^{i\cdot k\,2\pi/n}\right|
$$
The sum of a geometric series is recognized herein:
$$
\sum_{z\in S} z = \frac{1}{n} \frac{1-r^{n/2}}{1-r} \quad \mbox{with} \quad r = e^{i\cdot 2\pi/n}
$$
Hence:
$$
\left|\sum_{z\in S} z\right| = \frac{1}{n} \left|\frac{1-e^{i\cdot 2\pi/n\cdot n/2}}{1-e^{i\cdot 2\pi/n}}\right| =
\frac{1}{n} \left|\frac{2\cdot i}{e^{i\cdot \pi/n}\left(e^{-i\cdot \pi/n}-e^{i\cdot \pi/n}\right)\cdot i}\right| =
\frac{\pi/n}{\sin(\pi/n)}\frac{1}{\pi}
$$
And a well known limit tells us that
$$
\lim_{n\to\infty} \left|\sum_{z\in S} z\right| = \frac{1}{\pi}
$$
in concordance with the numerical experiments.
Don't get me wrong. The above is only a sketch of a proof. Quite some technicalities remain to be filled in.
The main part to be proved is: why should this very special case be relevant for the general case of arbitrary $\,z_k$ ?
(Though it's not uncommon that solutions of high symmetry are relevant for finding extreme values in a more general setting)