This question has been completely reformulated and a new property for the function $f_q$ has been added due to a series of helpful comments by fedja.
Consider the integral from quantum field theory due to F.A. Smirnov (see this MSE post): \begin{align} &\int_{-\infty}^{\infty}\prod_{j=1}^3\Gamma\Bigl(\frac 1 3 -\frac {\alpha-\beta_j}{2\pi i}\Bigr) \Gamma\Bigl(\frac 1 3 +\frac {\alpha-\beta_j}{2\pi i}\Bigr)(3\alpha-\sum\beta_m)e^{-\frac{\alpha}2}d\alpha\nonumber\\ &=\frac{(2\pi \Gamma(\frac 2 3))^2}{\Gamma(\frac 4 3)} \prod_{k\neq j}\Gamma\Bigl(\frac 2 3 -\frac {\beta_k-\beta_j}{2\pi i}\Bigr) e^{-\frac12\sum\beta_m}\sum e^{\beta_m},\quad |\text{Im}~\beta_j|<2\pi/3. \end{align} Some analytical and numerical calculations suggest that it has a $q$-analog of the form \begin{align} &\int_{-\infty}^\infty\prod_{j=1}^3\frac{\Gamma_q\left(\frac13+\frac{x-\beta_j}{2\pi i}\right)}{\Gamma_q\left(\frac23+\frac{x-\beta_j}{2\pi i}\right)}\frac{\sum q^{-\frac{\beta_m}{2 \pi i}}-\left(1+q^\frac13+q^{-\frac13}\right) q^{-\frac{x}{2 \pi i}}}{\prod_{m}\sin\left(\frac{\pi}3-\frac{x-\beta_m}{2 i}\right)}e^{-\frac{x}2}dx\\ &=\frac{-2\pi i}{(q;q)_{\infty }^9}\frac{\Gamma_q^2\left(\frac23\right)}{\Gamma_q\left(\frac43\right)}\frac{q^{5/9}}{(1-q)^2}e^{-\frac{1}{2} (\beta_1+\beta_2+\beta_3)} q^{-\frac{\beta_1+\beta_2+\beta_3}{6 \pi i}}\prod_{k\neq j}\Gamma_q\left(\frac23+\frac{\beta_j-\beta_k}{2\pi i}\right)\\ &\times\frac{\theta_q(q^{\frac{\beta_1-\beta_2}{2\pi i}})}{\sin\frac{\beta_1-\beta_2}{2 i}}\frac{\theta_q(q^{\frac{\beta_2-\beta_3}{2\pi i}})}{\sin\frac{\beta_2-\beta_3}{2 i}} \frac{\theta_q(q^{\frac{\beta_3-\beta_1}{2\pi i}})}{\sin\frac{\beta_3-\beta_1}{2 i}} \cdot f_{q_1}(\beta_1,\beta_2,\beta_3).\tag{5} \end{align} where $q_1=e^{-\frac{4 \pi ^2}{\ln(1/q)}}$, $\Gamma_q$ is the $q$-Gamma function and $$ \theta_q(z)=(q;q)_\infty(z;q)_\infty(q/z;q)_\infty. $$ The exact form of the function $f_q$ is unknown but there is evidence that it has a relatively simple closed form expression in terms of finite combination of theta functions.
Obviously $f_q(x,y,z)$ is symmetric in all three arguments $x,y,z$. It is known also that it has the following properties \begin{align} &f_q(x+c,y+c,z+c)=e^cf_q(x,y,z),\tag{1}\\ \\ &f_q\left(\frac{2\pi i}{3},0,z\right) \begin{aligned}[t] & =\left(e^z-1\right) \frac{\theta_q\left(e^{\frac{4 i \pi }{3}-z}\right)}{\theta_q\left(e^{-z}\right)}\\ & =\left(e^z-1\right)e^{-\frac{\pi i}{3}} q_1^{-\frac{1}{9}+\frac{z}{6 \pi i}}\frac{\theta_{q_1}\left(q_1^{\frac13+\frac{z}{2\pi i}}\right)}{\theta_{q_1}\left(q_1^{\frac{z}{2\pi i}}\right)} \end{aligned} \tag{2}\\ \\ &f_q(x,y,z)=e^x+e^y+e^z\\ &+\frac32\left(e^x+e^y+e^z\right)\left(1+e^{x-y}+e^{y-x}+e^{y-z}+e^{z-y}+e^{z-x}+e^{x-z}\right)q\\ &+i\frac{\sqrt{3}}{2}\left(e^{-x}+e^{-y}+e^{-z}\right)\left(e^{2x}+e^{2y}+e^{2z}\right)q+O(q^2) \tag{3} \\ &f_{q}(x,y,z)= e^{-\frac{\pi i}{3}} q_1^{-\frac29} \small{\frac{\left(e^z-e^y\right) q_1^{\frac{-2 x+y+z}{6 \pi i}}+\left(e^x-e^z\right) q_1^{\frac{x-2 y+z}{6 \pi i}}+\left(e^y-e^x\right) q_1^{\frac{x+y-2 z}{6 \pi i}}}{\Big(1-q_1^{\frac{x-y}{2 \pi i}}\Big) \Big(1-q_1^{\frac{z-x}{2 \pi i}}\Big) \Big(1-q_1^{\frac{y-z}{2 \pi i}}\Big)}\left(1+O(q_1^{\frac13})\right)}\tag{4} \\ \end{align}
Eq. $(2)$ is a condition that the difference of LHS and RHS of $(5)$ is an entire function. $(3)$ has been extracted from numerical calculations. $(4)$ has been found from exact $q$-series representation of the integral assuming $\beta_j\in\mathbb{R}$. However numerical calculations show that asymptotics $(4)$ holds if all differences $|\text{Im}~(\beta_j-\beta_k)|\leq \frac{2\pi}{3}$, and one can check that $(4)$ is consistent with $(2)$.
Q: Can anybody make a good guess from these data $1-4$ what this function $f_q(x,y,z)$ might be?
