以下は、crystal_MP_Ewald.py のソースコードを解析し、Sphinx（MyST）で問題なくビルド可能なMarkdownドキュメントです。

マデルングポテンシャル計算スクリプト crystal_MP_Ewald.py

概要

このスクリプト crystal_MP_Ewald.py は、エワルド法を用いて結晶のマデルングポテンシャルを計算するプログラムです。 指定された結晶の格子パラメータとサイト情報に基づき、実空間和、逆空間和、および自己項の3つの成分を合計して、最終的なポテンシャルを求めます。

機能

crystal_MP_Ewald.py は以下の主要な機能を提供します。

  • 初期化された格子パラメータとサイト情報（原子の種類、電荷、位置など）を使用します。

  • コマンドライン引数からエワルドパラメータ alpha と計算精度 prec を受け取ることができます。

  • 格子ベクトル、逆格子ベクトル、体積、および関連するメトリックテンソルを計算し表示します。

  • エワルドパラメータに基づき、実空間および逆空間の計算範囲（rdmax, G2max）を決定します。

  • 選択された中心サイトに対する実空間和 (UC1) を計算します。この項は実空間でのクーロン相互作用を表します。

  • 逆空間和 (UC2) を計算します。この項は逆格子空間でのクーロン相互作用を表し、高速フーリエ変換に似た形式です。

  • 自己相互作用項 (UC3) を計算します。これは原子自身の電場による自己エネルギーを補正する項です。

  • これら3つの項を合計し、マデルングポテンシャルおよびマデルング定数をJouleとeV単位で出力します。

  • 計算にかかった時間も合わせて表示されます。

非標準ライブラリ

このスクリプトは、以下の非標準ライブラリを使用しています。

  • numpy: 数値計算、特に配列操作や数学関数に利用されます。

  • matplotlib: グラフ描画ライブラリですが、このスクリプトでは matplotlib.pyplot がインポートされているものの、実際に描画処理は行われていません。

  • tkcrystalbase: 結晶学における基本的な計算（格子ベクトル、メトリックテンソル、体積、距離など）を提供するカスタムモジュールです。

使い方

プログラムはコマンドラインから以下の形式で実行します。エワルドパラメータ alpha と計算精度 prec はオプションで指定できます。

python crystal_MP_Ewald.py [alpha] [prec]

  • alpha: エワルドパラメータを指定します。省略された場合はコード内で定義されたデフォルト値 ew_alpha (0.3) が使用されます。

  • prec: 計算精度を指定します。省略された場合はコード内で定義されたデフォルト値 prec (1.0e-5) が使用されます。

例

デフォルト値で実行する場合：

python crystal_MP_Ewald.py

エワルドパラメータを 0.5、計算精度を 1.0e-6 で実行する場合：

python crystal_MP_Ewald.py 0.5 1.0e-6

入力

スクリプトは以下の情報を入力として使用します。

グローバル変数による初期設定

プログラムの冒頭で以下の変数が定義され、計算の初期設定として利用されます。

  • lattice_parameters: 格子定数と角度のリストです。

    • 例: [a, b, c, alpha, beta, gamma] (単位: オングストローム、度)

    • デフォルト値: [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0] (立方晶NaCl型構造に対応)

  • sites: 単位胞内の各サイトの情報を格納したリストです。各サイトは以下の要素を持つリストで表現されます。

    • ['原子名', 'サイトラベル', 原子番号, 原子質量, 電荷, 半径, '色', numpy.array([x, y, z])]

    • [x, y, z] は分数座標です。

    • デフォルト値: NaCl型構造の Na と Cl の8サイトが設定されています。

  • rmin: 同一サイトと判断する最小距離です。この距離よりも近いサイトは、実空間和の計算でスキップされます。

    • デフォルト値: 0.1

  • ew_alpha: エワルドパラメータです。

    • デフォルト値: 0.3

  • prec: 計算精度です (例: 1.0e-5)。

    • デフォルト値: 1.0e-5

コマンドライン引数

ew_alpha と prec は、コマンドライン引数によって上書きすることが可能です。詳細は「使い方」セクションを参照してください。

出力

スクリプトは計算結果を標準出力に表示します。出力される情報の主な項目は以下の通りです。

  • 格子パラメータと格子ベクトル:

    • 設定された格子パラメータ。

    • ax, ay, az の格子ベクトル（単位: A）。

    • メトリックテンソル gij（単位: A）。

    • 単位胞の体積（単位: A^3）。

  • 逆格子パラメータと逆格子ベクトル:

    • 逆格子パラメータ。

    • Rax, Ray, Raz の逆格子ベクトル（単位: A^-1）。

    • 逆格子メトリックテンソル Rgij（単位: A^-1）。

    • 逆格子単位胞の体積（単位: A^-3）。

  • エワルドパラメータ:

    • 使用された ew_alpha。

    • 設定された prec と norder。

  • 計算範囲:

    • 実空間和の最大距離 rdmax と、そのときの erfc(alpha*RDmax) の値。

    • 実空間和の繰返し回数の最大値 nrmax。

    • 実空間和で計算されるおおよそのサイト数 cal_N(real).

    • 逆空間和の最大 G^2 値 G2max と、そのときの exp(-pi^2 * G2max^2 / alpha^2) の値。

    • 逆空間和の繰返し回数の最大値 hgmax。

    • 逆空間和で計算されるおおよそのGベクトル数 cal_N(reciprocal).

  • 計算時間:

    • 実空間和の計算にかかった時間。

    • 逆空間和の計算にかかった時間。

    • 総計算時間。

  • マデルングポテンシャル:

    • 最終的なマデルングポテンシャル (MP) の値（Joule単位とeV単位）。

    • 実空間和 (UC1)、逆空間和 (UC2)、自己相互作用項 (UC3) の各成分がJouleとeV単位で示されます。

  • マデルング定数:

    • 計算されたマデルング定数。

コード解説

物理定数と数学定数

スクリプトの冒頭で、pi、pi2、torad、todeg、basee などの数学定数と、 h (プランク定数)、e (素電荷)、me (電子質量)、e0 (真空の誘電率) などの物理定数が定義されています。 特に、電荷の単位変換に使用される e2_4pie0 や、最終的なポテンシャル計算で使用される Ke (クーロン定数 e^2 / (4 * pi * e0)) が定義されています。

初期データ設定

lattice_parameters と sites グローバル変数に、計算対象となる結晶の初期データが設定されています。 これは、NaCl型構造に対応する立方晶の格子とサイト配置の例です。

usage() 関数

def usage():
    """
    概要:
        プログラムの正しい使用方法を標準出力に表示します。
    詳細説明:
        プログラムをコマンドラインから実行する際の引数のフォーマットと例を示します。
    """
    print("")
    print("Usage: python {} alpha prec".format(argv[0]))
    print("   ex: python {} {} {}".format(argv[0], ew_alpha, prec))
    print("")

この関数は、スクリプトの正しい使用方法とコマンドライン引数の例を標準出力に表示します。

terminate() 関数

def terminate():
    """
    概要:
        使用方法を表示した後、プログラムを終了します。
    詳細説明:
        usage() 関数を呼び出し、その後にプログラムを強制終了します。
    """
    usage()
    exit()

usage() 関数を呼び出した後、プログラムを終了させるための補助関数です。

main() 関数

def main():
    """
    概要:
        エワルド法によりマデルングポテンシャルを計算し、結果を表示します。
    詳細説明:
        この関数は、設定された格子パラメータとサイト情報に基づき、
        エワルド法を用いて結晶のマデルングポテンシャルを計算し、その結果を標準出力に表示します。
        具体的な計算手順は以下の通りです。
        ... (詳細説明は上記の「機能」と「出力」セクションにまとめられています)
    """
    # ...

プログラムの主要なロジックを実装している関数です。以下の手順でマデルングポテンシャルの計算と結果の表示を行います。

  1. 格子情報と逆格子情報の計算と表示: cal_lattice_vectors()、cal_metrics()、cal_volume() (tkcrystalbase モジュールから) を使用して、格子ベクトル、メトリックテンソル、体積を計算し表示します。 同様に、cal_reciprocal_lattice_vectors()、cal_reciprocal_lattice_parameters()、cal_metrics() (tkcrystalbase モジュールから) を使用して、逆格子関連の情報を計算し表示します。

  2. エワルドパラメータと計算範囲の決定: コマンドライン引数またはデフォルトで設定された ew_alpha と prec を表示します。 これらの値に基づき、実空間和の最大距離 rdmax と逆空間和の最大Gベクトル二乗値 G2max を計算します。 さらに、それぞれの計算範囲に対応する最大繰返し回数 nrmax および hgmax を推定し表示します。 これらの計算には erfc() や log10()、log() などの数学関数が使用されます。

  3. 実空間和 (UC1) の計算: UC1 = 0.0 3重の for ループ (ix, iy, iz) を用いて、中心サイト sites[0] の周りの周期的な単位胞を探索します。 各単位胞内の全サイト j との距離 rij を distance() 関数 (tkcrystalbase モジュールから) で計算します。 もし rij < rmin であれば、同一サイトと判断して計算をスキップします。 erfc(ew_alpha * rij) / (rij * 1.0e-10) を各サイトの電荷 qj と乗算し、UC1 に加算します。この結果はeV単位に変換されます。

  4. 逆空間和 (UC2) の計算: UC2 = 0.0 3重の for ループ (h, k, l) を用いて、逆格子空間のGベクトルを探索します。 G^2 値を distance2() 関数 (tkcrystalbase モジュールから) で計算します。 もし G^2 == 0.0 (中心Gベクトル) または G^2 > G2max であれば、計算をスキップします。 各サイト j の電荷 qj と位置 pos_j を用いて構造因子に似た和 cossum_j と sinsum_j を計算します。 中心サイト i の位相 phi_i とこれらの和から fcal を計算し、exp(-Kexp * G2) / (G2 * 1.0e+20) と Krec を乗算して UC2 に加算します。

  5. 自己相互作用項 (UC3) の計算: UC3 = qi * 2.0 * (ew_alpha * 1.0e10) / sqrt(pi) 中心サイト i の電荷 qi とエワルドパラメータ ew_alpha を用いて計算されます。これはエワルド法で導入される自己相互作用を打ち消すための補正項です。

  6. 結果の集計と表示: MP = UC1 + UC2 - UC3 として最終的なマデルングポテンシャルを計算します。 各計算フェーズにかかった時間とともに、MP をJoule単位とeV単位で表示します。 最後に、マデルングポテンシャルからマデルング定数を計算し表示します。

メインガード

if __name__ == '__main__':
    argv = sys.argv
    narg = len(argv)
    if narg >= 2:
        ew_alpha = float(argv[1])
    if narg >= 3:
        prec = float(argv[2])

    main()

スクリプトが直接実行された際に、sys.argv を解析してコマンドライン引数から ew_alpha と prec を取得し、main() 関数を呼び出します。

注意事項

コードからは、特定の利用上の制限や特別な注意事項は確認できません。一般的なPythonスクリプトの実行環境があれば動作すると考えられます。