In MNDO and its descendants (Sewart's PMx models, Austin's AM1, MNDO/d, etc.), the nonvanishing two-center two-electron "coulomb integrals" are evaluated via a classical multipole expansion in a local coordinate system, and then transformed into global coordinates via "standard rotational transformations". How exactly are these transformations carried out?
A similar question would be the methodology behind the resonance integrals in the orthogonalised methods OMx, which are similarly evaluated in a diatomic coordinate system and then transformed to the global coordinate system.