Densities

In order to calculate the potentials in Eq. (46) we have to determine all secondary densities:

$$\rho_{mI,nLvJ,J'M'}(\vec{r}) = [D_{mI}\rho_{nLvJ}(\vec{r})]_{J'M'} = [D_{mI}[K_{nL}\rho_v(\vec{r}_1,\vec{r}_2)]_J]_{J'M'},$$ (67)

where $D_{mI}$ is built by coupling gradients $\nabla_1+\nabla_2$ and $K_{nL}$ is built by coupling gradients $(\nabla_1-\nabla_2)/2i$. Here and below we understand that $\vec{r}=\vec{r}_1=\vec{r}_2$ is set after performing all the differentiations.

In analogy with Eq. (55), we can split these operators as

$$K_{nLM_L} \!\! = \sum_{rRr'R'} K^{nL}_{rRr'R'}\sum_{M_RM'_R} C^{LM_L}_{RM_RR'M'_R}D^{(1)}_{rRM_R}D^{(2)}_{r'R'M'_R},$$ (68)

$$D_{mIM_I} \!\! = (2i)^m\sum_{pPp'P'} (-1)^{p'} K^{mI}_{pPp'P'}\sum_{M_PM'_P} C^{IM_I}_{PM_PP'M'_P}D^{(1)}_{pPM_P}D^{(2)}_{p'P'M'_P},$$ (69)

where $n=r+r'$ and $m=p+p'$. Then the density (68) reads

$$\rho_{mI,nLvJ,J'M'}(\vec{r}) \!\!\!\! = \sum_{M_IM_J}C^{J'M'}_{IM_IJM_J}\sum_{M_L\mu}C^{JM_J}_{LM_Lv\mu}$$ (70)

$$\times(2i)^m\sum_{pPp'P'} (-1)^{p'} K^{mI}_{pPp'P'}\sum_{M_PM'_P} C^{IM_I}_{PM_PP'M'_P}D^{(1)}_{pPM_P}D^{(2)}_{p'P'M'_P}$$

$$\times \sum_{rRr'R'} K^{nL}_{rRr'R'}\sum_{M_RM'_R} C^{LM_L}_{RM_RR'M'_R}D^{(1)}_{rRM_R}D^{(2)}_{r'R'M'_R} \rho_{v\mu}(\vec{r}_1,\vec{r}_2).$$

We now introduce coefficients $D^{uU}_{rRpP}$, which allow for expressing products of derivatives as:

$$D^{(1)}_{rRM_R}D^{(1)}_{pPM_P} = \sum_{uUM_U} C^{UM_U}_{RM_RPM_P} D^{uU}_{rRpP} D^{(1)}_{uUM_U} ,$$ (71)

$$D^{(2)}_{r'R'M'_R}D^{(2)}_{p'P'M'_P} = \sum_{u'U'M'_U} C^{U'M'_U}_{R'M'_RP'M'_P} D^{u'U'}_{r'R'p'P'} D^{(2)}_{u'U'M'_U} ,$$ (72)

where $u=r+p$ and $u'=r'+p'$, that is, $u+u'=m+n$. These coefficients can be calculated by using methods outlined in A. At N$^3$LO, only 91 coefficients $D^{uU}_{rRpP}$ are needed, so they can easily be precalculated and stored. The sum of products of four Clebsh-Gordan coefficients can now be recoupled (see Eq. 8.7(20) in Ref. [6]) as:

$$\sum_{M_PM'_P}\sum_{M_RM'_R} C^{IM_I}_{PM_PP'M'_P} C^{LM_L}_{RM_RR'M'_R} C^{UM_U}_{RM_RPM_P} C^{U'M'_U}_{R'M'_RP'M'_P}$$

$$= \sqrt{(2I+1)(2L+1)(2U+1)(2U'+1)}$$

$$\times \sum_{WM_W} \left\{\begin{array}{rrr} P' & P & I \\ R' ... ...L \\ U' & U & W \end{array}\right\} C^{WM_W}_{UM_UU'M'_U}C^{WM_W}_{LM_LIM_I}$$ (73)

Subsequently, the sum of products of three Clebsh-Gordan coefficients can be recoupled (see Eq. 8.7(12) in Ref. [6]) as:

$$\sum_{M_JM_LM_I} C^{WM_W}_{LM_LIM_I} C^{J'M'}_{IM_IJM_J} C^{JM_J}_{LM_Lv\mu}$$

$$= (-1)^{W+v-J'}\sqrt{(2W+1)(2J+1)} \left\{\begin{array}{rrr} L & I & W \\ J' & v & J \end{array}\right\} C^{J'M'}_{WM_Wv\mu} .$$ (74)

The last two remaining Clebsh-Gordan coefficients can be absorbed in the following definition of the coupled derivative of the density:

$$\rho^{uUu'U'W}_{vJ'M'}(\vec{r}) = [[D^{(1)}_{uU}D^{(2)}_{u'U'}]_W\rho_v(\vec{r}_1,\vec{r}_2)]_{J'M'}$$ (75)

$$= \sum_{M_W\mu} C^{J'M'}_{WM_Wv\mu} \sum_{M_UM'_U}C^{WM_W}_{UM_UU'M'_U} D^{(1)}_{uUM_U}D^{(2)}_{u'U'M'_U} \rho_{v\mu}(\vec{r}_1,\vec{r}_2) ,$$

which finally gives

$$\rho_{mI,nLvJ,J'M'}(\vec{r}) \!\! = \sum_{uUu'U'W} A^{uUu'U',W}_{mI,nLvJ,J'} \rho^{uUu'U'W}_{vJ'M'}(\vec{r}) ,$$ (76)

where coefficients $A^{uUu'U',W}_{mI,nLvJ,J'}$ result from summing up all intrinsic indices:

$$A^{uUu'U',W}_{mI,nLvJ,J'} \!\! = (2i)^m\sum_{pPp'P'} \sum_{rRr'R'} (-1)^{p'} K^{mI}_{pPp'P'} K^{nL}_{rRr'R'} D^{uU}_{rRpP}D^{u'U'}_{r'R'p'P'}$$

$$\times \sqrt{(2I+1)(2L+1)(2U+1)(2U'+1)} \left\{\begin{array}{rrr} P' & P & I \\ R' & R & L \\ U' & U & W \end{array}\right\}$$

$$\times(-1)^{W+v-J'}\sqrt{(2W+1)(2J+1)} \left\{\begin{array}{rrr} L & I & W \\ J' & v & J \end{array}\right\} .$$ (77)

At N$^3$LO, only 3138 coefficients $A^{uUu'U',W}_{mI,nLvJ,J'}$ are needed, so they can easily be precalculated and stored.