Firstly, I want to show that, in the optimal configuration, the points must lie on the boundary. Let $z_1,\ldots, z_n\in\mathbb{C}$ be the complex numbers that represent the $n$ points in this optimal configuration and suppose, without loss of generality, that $z_n\notin\partial D$. Consider the function $f: D\rightarrow\mathbb{C}$, $f(z)=\prod_{i<n}(z-z_i)$. Since $f$ is holomorphic and non-constant, by the maximum modulus principle, $|f(z)|$ must attain its maximum on the boundary $\partial D$. It follows that there exist $z'\in\partial D$ : $|f(z')|>|f(z_n)|$, which contradicts the optimality of the configuration. Therefore, assume $z_1,\ldots, z_n$ lie on the boundary of the unit disk and assume also that $z_i=\cos\theta_i+i\sin\theta_i$, with $0\le\theta_1<\ldots<\theta_n<2\pi$. Let $\phi_i=\theta_{i+1}-\theta_i$ for $i=1,\ldots, n$ (with $\theta_{n+1}\equiv \theta_1$). It follows that $\sum_i\phi_i=2\pi$. Then $|z_{i+1}-z_i|=\sqrt{2-2\cos\phi_i}=2\sin\frac{\phi_i}{2}$, and the following holds: \[\prod_i|z_{i+1}-z_i|=2^n\prod_i\sin\frac{\phi_i}{2}\le 2^n\left(\sum_i \frac{1}{n}\sin\frac{\phi_i}{2}\right)^n\le 2^n \sin\left(\frac{1}{2n}\sum_i\phi_i\right)^n=\] \[=\left(2\sin\frac{\pi}{n}\right)^n.\] where in the first inequality I used AM-GM, while the second is Jensen's applied to the sine function, which is concave in $[0,\pi]$. In the same way, we can show that \[\prod_i|z_{i+2}-z_i|=2^n\prod_i\sin\frac{\phi_i+\phi_{i+1}}{2}\le \left(2\sin\frac{2\pi}{n}\right)^n,\] and it's not difficult to guess the general inequality, for $k\le \left\lfloor\frac{n}{2}\right\rfloor$, \[\prod_i|z_{i+k}-z_i|\le \left(2\sin\frac{k\pi}{n}\right)^n.\] From this follows that, \[\prod_{k=1}^{ \left\lfloor\frac{n}{2}\right\rfloor}\prod_i|z_{i+k}-z_i|\le\prod_{k=1}^{ \left\lfloor\frac{n}{2}\right\rfloor} \left(2\sin\frac{k\pi}{n}\right)^n.\] The optimal configuration is clearly reached when the points are placed on the vertices of a regular $n$-gon.

Here's another solution. The answer is $n^{n / 2}$, which may be achieved by letting the points be the vertices of an inscribed regular $n$-gon. Now, we show that the product $P$ cannot exceed $n^{n / 2}$. Let the points have complex coordinates $a_1, \dots, a_n$, with $|a_k| \le 1$ for all $k$. Then, by the Vandermonde determinant \[P = \prod_{i < j} |a_i - a_j| = |\det A|,\]where \[A = \begin{bmatrix} 1 & a_1 & \dots & a_1^{n - 1} \\ 1 & a_2 & \dots & a_2^{n - 1} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & a_n & \dots & a_n^{n - 1} \end{bmatrix}.\]In particular, note that each entry of $A$ has magnitude at most $1$. We will show that $|\det A| \le n^{n / 2}$ for any $n \times n$ matrix $A$ with entries in the unit disk. In fact, if the column vectors comprising $A$ are $v_1, \dots, v_n$, we show that $|\det A| \le |v_1| \cdots |v_n|$, which implies the conclusion. By scaling each column, we may suppose that $|v_1| = \dots = |v_n| = 1$, in which case we wish to show $|\det A| \le 1$. Let $M = A^{\dagger}A$, where $A^{\dagger}$ is the conjugate transpose. It follows that $M$ is Hermitian and positive semidefinite, as $v^{\dagger} A^{\dagger} A v = x^{\dagger} x \ge 0$. Thus the eigenvalues $\lambda_1, \dots, \lambda_n$ of $M$ are nonnegative. Now, since $\det M = \lambda_1 \cdots \lambda_n$ and $\operatorname{tr} M = \lambda_1 + \dots + \lambda_n = \overline{v_1} \cdot v_1 + \dots + \overline{v_n} \cdot v_n = n$, we deduce \[|\det A|^2 = |\det M| = |\lambda_1 \cdots \lambda_n| \le \left | \frac{\lambda_1 + \dots + \lambda_n}{n} \right |^n = 1\]by AM-GM, and thus $|\det A| \le 1$ as desired.