Is there a way to obtain solution-phase dielectric constants?

It's possible to estimate solution-phase dielectric constant from a molecular dynamics simulation using this formula:

$$epsilon_{r} = 1 + frac{4pi}{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.