crystal_MP_simple.py ドキュメント
.. toctree:: :maxdepth: 2 :hidden:
self
概要
crystal_MP_simple.py は、結晶のマデルングポテンシャルを単純総和法で計算し、結果をプロットするPythonスクリプトです。
スクリプトの説明
このスクリプトは、指定された結晶構造（格子定数とサイト情報）に基づき、単純総和法を用いてマデルングポテンシャルを計算します。計算されたポテンシャルは、中心イオンからの距離の関数としてプロットされます。
tkcrystalbase.py モジュールに依存します。
関連リンク: crystal_MP_simple_usage
非標準ライブラリ
このスクリプトは以下の非標準ライブラリを使用します。
numpy(np): 数値計算のためのライブラリ。numpy.linalg(la): 線形代数演算のためのモジュール。matplotlib.pyplot(plt): グラフ描画のためのライブラリ。tkcrystalbase: 結晶構造に関するユーティリティ関数を提供するカスタムモジュール。
定数
スクリプト内で定義されている主要な定数とその値、単位は以下の通りです。
pi: 3.14159265358979323846pi2:piの2倍torad: 0.01745329251944 (rad/deg)todeg: 57.29577951472 (deg/rad)basee: 2.71828183 (自然対数の底)h: 6.6260755e-34 (Js, プランク定数)h_bar(hbar): 1.05459e-34 (Js, ディラック定数)c: 2.99792458e8 (m/s, 光速)e: 1.60218e-19 (C, 素電荷)me: 9.1093897e-31 (kg, 電子の質量)mp: 1.6726231e-27 (kg, プロトンの質量)mn: 1.67495e-27 (kg, 中性子の質量)u0: 4.0 * 3.14 * 1e-7 (Ns^2 C^-2, 真空の透磁率)e0: 8.854418782e-12 (C^2 N^-1 m^-2, 真空の誘電率)e2_4pie0: 2.30711e-28 (Nm^2,e^2 / (4 * pi * e0))a0: 5.29177e-11 (m, ボーア半径)kB: 1.380658e-23 (JK^-1, ボルツマン定数)NA: 6.0221367e23 (mol^-1, アボガドロ定数)R: 8.31451 (J/K/mol, 理想気体定数)F: 96485.3 (C/mol, ファラデー定数)g: 9.81 (m/s^2, 重力加速度)
グローバル変数
スクリプト内で定義されているグローバル変数は以下の通りです。
lattice_parameters: 結晶の格子定数を格納するリスト。[a, b, c, alpha, beta, gamma]の形式で、単位はオングストローム (angstrom) と度 (degree) です。 デフォルト値は[5.62, 5.62, 5.62, 90.0, 90.0, 90.0]（立方晶）です。コメントアウトされた行には[5.62, 5.62, 5.62, 60.0, 60.0, 60.0]（菱面体晶）の例もあります。sites: 結晶サイト情報を格納するリストです。各要素は以下の情報を持つリストです。原子名 (例:
'Na','Cl')サイトラベル (例:
'Na1','Cl1')原子番号
原子量
電荷 (例:
+1.0,-1.0)半径
色 (プロット用)
位置 (
numpy.ndarray型の分数座標)
rmin: マデルングポテンシャル計算およびプロットの開始距離 (angstrom)。デフォルト値は0.1です。rmax: マデルングポテンシャル計算およびプロットの最大距離 (angstrom)。デフォルト値は100.0です。コマンドライン引数で上書き可能です。nr:rminからrmaxまでのプロット点数。デフォルト値は101です。コマンドライン引数で上書き可能です。figsize:matplotlibでの図のサイズ (幅, 高さ) を指定するタプル。デフォルト値は(6, 6)です。rstep:rminとrmaxの間の距離ステップ幅。(rmax - rmin) / (nr - 1)で計算されます。
コマンドライン引数
このスクリプトは以下のコマンドライン引数を受け取ります。
rmax(第1引数): マデルングポテンシャルを計算する最大距離を浮動小数点数で指定します。指定しない場合、グローバル変数rmaxのデフォルト値 (100.0) が使用されます。nr(第2引数):rminからrmaxまでの計算点数を整数で指定します。指定しない場合、グローバル変数nrのデフォルト値 (101) が使用されます。
入出力仕様
入力ファイル: このスクリプトは、外部の入力ファイルを読み込みません。結晶構造情報（
lattice_parametersとsites）はスクリプト内部で直接定義されています。出力:
標準出力に、格子定数、格子ベクトル、メートルテンソル、単位格子体積、逆格子ベクトル、逆格子メートルテンソル、逆格子体積、計算範囲 (
nmax) などの情報が出力されます。距離
rとそれに対応するマデルングポテンシャル (MP) の値が表形式で標準出力されます。matplotlibを用いて、距離r(オングストローム) に対する静電ポテンシャル (電子ボルト) のグラフが画面に表示されます。ファイルへの保存はコードからは確認できません。
主要な関数
usage()
概要: スクリプトの正しい使用方法を表示します。
詳細説明: コマンドライン引数 rmax と nr の指定方法をユーザーに示します。この関数内で argv[0] が参照されていますが、argv は if __name__ == '__main__': ブロック内のローカル変数であり、グローバル変数として定義されていません。このため、usage() が直接または terminate() を介して呼び出された場合、NameError が発生する可能性があります。
引数: :returns: なし :rtype: None
terminate()
概要: エラー発生時に使用方法を表示してスクリプトを終了します。
詳細説明: usage() 関数を呼び出して使用方法を出力した後、システムを終了します。
引数: :returns: なし (スクリプトが終了するため) :rtype: None
draw_box(ax, aij, nrange, color='black')
概要: 結晶の単位格子境界を3Dプロットに描画します。
詳細説明: 単位格子の辺を黒線（または指定された色）で描画します。この関数は draw_unitcell() から呼び出されますが、nrange は現在の実装では使用されていません。
引数:
:param ax: matplotlib の3D軸オブジェクト。
:type ax: matplotlib.axes._subplots.Axes3DSubplot
:param aij: (3, 3)の ndarray 、格子ベクトル a, b, c を表す。
:type aij: numpy.ndarray
:param nrange: 描画する単位格子の範囲。 [[xmin, xmax], [ymin, ymax], [zmin, zmax]] の形式。（現状未使用）
:type nrange: list of list of float
:param color: 描画する線の色。デフォルトは 'black'。
:type color: str
戻り値:
:returns: なし
:rtype: None
draw_unitcell(ax, sites, aij, nrange, color='black')
概要: 結晶の単位格子とその中の原子を3Dプロットに描画します。
詳細説明: draw_box() 関数を呼び出して単位格子を描画し、その後、sites リスト内の原子を分数座標からデカルト座標に変換してプロットします。nrange は描画する単位格子の範囲を指定しますが、このスクリプトの main() 関数では現在呼び出されていません。原子の描画サイズは kr * r で計算されますが、変数 kr はこのスクリプトのコード内では定義されておらず、tkcrystalbase モジュールから提供されている可能性が高いですが、コードからはその定義を確認できません。
引数:
:param ax: matplotlib の3D軸オブジェクト。
:type ax: matplotlib.axes._subplots.Axes3DSubplot
:param sites: サイト情報のリスト。各サイトは [atom_name, site_label, atomic_number, atomic_mass, charge, radius, color, position] の形式。
:type sites: list of list
:param aij: (3, 3)の ndarray 、格子ベクトル a, b, c を表す。
:type aij: numpy.ndarray
:param nrange: 描画する単位格子の範囲。 [[xmin, xmax], [ymin, ymax], [zmin, zmax]] の形式。
:type nrange: list of list of int
:param color: 単位格子を描画する線の色。デフォルトは 'black'。
:type color: str
戻り値:
:returns: なし
:rtype: None
main()
概要: マデルングポテンシャルの計算と結果のプロットを実行します。
詳細説明:
格子定数 (lattice_parameters) から tkcrystalbase モジュール内の関数 (cal_lattice_vectors(), cal_metrics(), cal_volume(), cal_reciprocal_lattice_vectors(), cal_reciprocal_lattice_parameters()) を用いて、格子ベクトル、メートルテンソル、単位格子体積、逆格子ベクトル、逆格子メートルテンソル、逆格子体積を計算し、その情報を表示します。
次に、rmax と nr に基づいて、原点にあるイオン (sites[0]) に対するマデルングポテンシャルを単純総和法で計算します。計算範囲は rmax から導出される nxmax, nymax, nzmax の範囲で格子点を走査します。各イオン間の距離は tkcrystalbase.distance() 関数で計算されます。ポテンシャルは Ke * q1 / (r * 1.0e-10) / e の式で計算され、結果は eV 単位に変換されます。
計算されたポテンシャルの微分値 (MPdiff) が距離 r の各ビンに累積され、その後、逐次的に合計されてマデルングポテンシャル (MP) の最終値が求められます。
計算されたポテンシャルは、距離 r に対するグラフとして matplotlib で表示されます。x 軸は 'r / angstrom'、y 軸は 'Electrostatic potential / eV' となります。
最後に、terminate() 関数が呼び出され、スクリプトが終了します。
戻り値: :returns: なし :rtype: None
使用例
crystal_MP_simple.py を実行する際のコマンドライン引数の例を以下に示します。
python crystal_MP_simple.py
この場合、rmax はデフォルトの 100.0、nr はデフォルトの 101 が使用されます。
rmax と nr を指定する場合:
python crystal_MP_simple.py 50.0 201
この例では、最大距離を 50.0 オングストローム、点数を 201 に設定して計算を実行します。