It's possible to estimate solution-phase dielectric constant from a molecular dynamics simulation using this formula:
$$\epsilon_{r} = 1 + \frac{4\pi}{3Vk_{B}T}(\langle \mathbf{P}^{2} \rangle - \langle \mathbf{P} \rangle^{2})$$
Where $V$ is the volume, $k_{B}$ is Boltzmann's constant, $T$ is temperature, and $P$ is the dipole moment defined as: $\mathbf{P} = \sum_{i} \vec{\mu}_{i}$ the summation of molecular dipole moments.
In the absence of any external electric field (which I assume is the case here), from electrostatics, you have:
$$\mathbf{P}(\mathbf{r}) = \chi \int_{\Omega} \mathbf{T}(\mathbf{r}-\mathbf{r}^{'})\cdot \mathbf{P}(\mathbf{r}^{'}) d^{3} \mathbf{r}^{'}$$
$\chi$ is the susceptibility, which is unknown here. Also, $\mathbf{T}$ is the dipole-dipole tensor defined as:
$$T_{ij} =
\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}(-\ln(r))$$
Now if you replace the integral with a summation and replace $\mathbf{P}(\mathbf{r})$ with the discretized dipole moment at molecular locations shown as $\mathsf{P}$ and the discretized dipole-dipole tensor (matrix $\mathsf{T}$), you have:
$$\mathsf{P} = \chi \mathsf{T} \cdot \mathsf{P}$$
or:
$$\mathsf{T} \cdot \mathsf{P} = \frac{1}{\chi} \mathsf{P}$$
This is an eigenvalue problem where you know dipole-dipole tensor $\mathsf{T}$, while the eigenvector (dipole moment $\mathsf{P}$) and eigenvalue (susceptibility $\chi$) are unknowns. You can solve this eigenvalue problem for your system and then you'll get the dielectric constant by estimating the fluctuation of dipole moment.