molecules. The Hamiltonian of the system, $$\begin{aligned}
H(\vec{r}_1,\vec{r}_2,\vec{R}) = T(\vec{r}_1) + T(\vec{r}_2) + T(\vec{R}) + V(\vec{r}_1,\vec{r}_2,\vec{R}) ,
\label{HAM}\end{aligned}$$ is composed of a radial kinetic energy term $T(\vec{R})$ describing the center-of-mass motion, two kinetic energy terms $T(\vec{r}_1)$ and $T(\vec{r}_2)$ for each diatomic molecule and the PES function, $$\begin{aligned}
V(\vec{r}_1,\vec{r}_2,\vec{R}) = U(\vec{r}_1,\vec{r}_2,\vec{R}) + V(\vec{r}_1) + V(\vec{r}_2) ,
\label{Potential}\end{aligned}$$ describing the interactions between the four atoms. The terms $V(\vec{r}_1)$ and $V(\vec{r}_2)$ in Eq.  are the interaction energy potentials of the two H$_2$ molecules and $U(\vec{r}_1,\vec{r}_2,\vec{R})$ is the interaction energy potential between the two molecules, which vanishes at large molecule - molecule separations, $$\begin{aligned}
\lim_{R\to\infty}^{} U(\vec{r}_1,\vec{r}_2,\vec{R}) \to 0. \end{aligned}$$ The angular dependence of the interaction potential may be expanded as $$\begin{aligned}
U(\vec{r}_1,\vec{r}_2,\vec{R}) = \sum^{}_{\lambda} \ 
A_{\lambda}(r_1,r_2,R)  Y_{\lambda}(\hat{r}_1,\hat{r}_2,\hat{R}) ,
\label{POTLAMBDA}\end{aligned}$$ with $$\begin{gathered}
Y_{\lambda}(\hat{r}_1,\hat{r}_2,\hat{R}) = 
\sum^{}_{m_{\lambda}} \big\langle  \lambda_1 m_{\lambda_1} \lambda_2 m_{\lambda_2} \big| \lambda_{12} m_{\lambda_{12}} \big\rangle \\
Y_{\lambda_1 m_{\lambda_1}}(\hat{r}_1) Y_{\lambda_2 m_{\lambda_2}}(\hat{r}_2) Y^*_{\lambda_{12} m_{\lambda_{12}}}(\hat{R}) ,\end{gathered}$$ where $\lambda \equiv \lambda_1 \lambda_2 \lambda_{12}$ and $m_{\lambda} \equiv m_{\lambda_1} m_{\lambda_2} m_{\lambda_{12}}$. The notation $\big\langle  \lambda_1 m_{\lambda_1} \lambda_2 m_{\lambda_2} \big| \lambda_{12} m_{\lambda_{12}} \big\rangle$ represents a Clebsch–Gordan coefficient. Because the Hamiltonian defined in Eq.  is invariant under a rotation in space, the total angular momentum $\vec J$ and its projection $J_z$ on a space-fixed 