I am interested in calculating the transition dipole moment (TDM) from the information from two wavefunctions of different states. This is somewhat similar to calculating the molecular dipole moment which was previously answered here:
How to calculate molecular dipole moment from a known wavefunction?
The information I have is the molecular orbital coefficients of alpha and beta from states 1 and 2 and the multipole matrix. In considering just one Cartesian direction X, calculating the dipole moment along X would be something like
$$\mathbf{P_a} = \mathbf{C_a} \cdot \mathbf{C_a^T}$$ $$\mathbf{P_b} = \mathbf{C_b} \cdot \mathbf{C_b^T}$$ $$\mathbf{P} = \mathbf{P_a} + \mathbf{P_b}$$ $$\mathbf{\mu_{matrix}} = \mathbf{P} \cdot \mathbf{\text{multipoleMatrixX}}$$ $$\mu_{electronic} = trace(\mathbf{\mu_{matrix}})$$
$$\mu_{nuclear} = \sum_{\text{all atoms}}(Z_{nucleus} * \overrightarrow{r_x})$$
$$\mu_{X} = - \mu_{electronic} + \mu_{nuclear}$$
where $\mathbf{C_a}$ and $\mathbf{C_b}$ are the occupied alpha and beta molecular orbital coefficients, $\mathbf{P_a}$ and $\mathbf{P_b}$ are the alpha and beta density, $\mathbf{P}$ is the total density, and $\mathbf{\text{multipoleMatrixX}}$ is the multipole matrix for the X direction.
I might think that to calculate the transition dipole moment I might just change the densities to the transition densities doing something like I have show below. The transition denisty is based on page 10 of the article here:
$$\mathbf{S_a} = \mathbf{C_{a1}} \cdot \mathbf{C_{a2}^T}$$ $$\mathbf{S_b} = \mathbf{C_{b1}} \cdot \mathbf{C_{b2}^T}$$ $$\mathbf{P_a} = \mathbf{C_{a1}} \cdot \mathbf{C_{a2}^T} \cdot \mathbf{S_a^{-1}}$$ $$\mathbf{P_b} = \mathbf{C_{b1}} \cdot \mathbf{C_{b2}^T} \cdot \mathbf{S_b^{-1}}$$ $$\mathbf{P} = \mathbf{P_a} + \mathbf{P_b}$$ $$\mathbf{\mu_{matrix}} = \mathbf{P} \cdot \mathbf{\text{multipoleMatrixX}}$$ $$\mu_{electronic} = trace(\mathbf{\mu_{matrix}})$$ $$\mu_{nuclear} = \sum_{\text{all atoms}}(Z_{nucleus} * \overrightarrow{r_x})$$
$$\mu_{X,TDM} = - \mu_{electronic} + \mu_{nuclear}$$
where $\mathbf{C_{a1}}$ and $\mathbf{C_{a2}}$ are the occupied alpha molecular orbital coefficients from state 1 and 2, $\mathbf{C_{b1}}$ and $\mathbf{C_{b2}}$ are the occupied beta molecular orbital coefficients from state 1 and 2, $\mathbf{S_a}$ and $\mathbf{S_b}$ are the alpha and beta overlap matricies, $\mathbf{S_a^{-1}}$ and $\mathbf{S_b^{-1}}$ are the inverse of the overlap matrices.
I am unsure how to calculate the TDM since it is in disagreement with known calculated TDM values; some of these transitions should have a TDM of zero due to symmetry considerations. I have been unable to find any discussion on calculating the TDM from two wavefunctions, and most discussions I have seen are in the context of some perturbation on the ground state (i. e. TD or CI). I would be grateful if anyone could provide some suggestions or references to look at.