TB法でエネルギーバンド計算

以前p系軌道の解析積分の計算結果を載せていたのですが、
あれは原子用に原点中心にした簡易バージョンです。
分子や結晶のエネルギーを計算する場合は、多体中心積分をやらないと駄目です。
s軌道系なら結構簡単に計算できますが、p軌道系以上になると途端に複雑になります。
こういったガウス軌道の積分方法は文献[1]が詳しいです。ただ、p軌道以上の
積分結果は「とても複雑になる」としか書かれていないです。
これといった簡単に計算できるアルゴリズムもないのかな?

ただし、計算の指針は書かれいます。
==============================
sガウス軌道関数を
g_{s}(r)=\exp\left(-\alpha\left(r-R\right)^{2}\right)
とする。R_{x}に関して微分をするとp系の軌道関数が得られる。
g_{px}(r)=xg_{s}(r)=\displaystyle \frac{1}{2\alpha}\frac{\partial g_{s}(r)}{\partial R_{x}}+R_{x}g_{s}(r)
となるので、各種解析積分は例えば以下のようにs系解析積分微分して計算できる。
\displaystyle \left\langle px,\alpha,R_{A}|T|px,\beta,R_{B}\right\rangle=\frac{1}{4\alpha\beta}\frac{\partial}{\partial R_{Ax}}\frac{\partial}{\partial R_{Bx}}\left\langle s|T|s\right\rangle+\cdots
==============================
というものです。確かにこれで機械的に計算できそうです。
ただ、2電子積分となると4階の微分操作が必要で、やっぱり面倒&複雑になりそうです。
こういうのはmathematicaで結果を展開するといいかもしれません。
一応、手動でp軌道の2中心のクーロン積分\left\langle px,\alpha,R_{A}|\frac{1}{|r-R_{C}|}|px,\beta,R_{B}\right\rangle)を計算しましたが、これで力つきました。

しかし、頑張って厳密にp軌道系を計算するようなことをせずとも、
実際にはp軌道系をs軌道系の和で近似する方法がとられているようです。

おまけ:
丸善分子キットで単純立方格子を作ってみた。
微妙に歪んでいるのは買ったのが有機系キット(無機系は高い)らしく
うまく組み立てられず、かなり無理して作ってます。

参考:
[1]新しい量子化学―電子構造の理論入門〈上〉
[2]金属電子論〈上〉
[3]コンピュータでみる固体の中の電子―バンド計算の基礎と応用