以下は、H1s_HF_LDA.py プログラムの解析に基づいたMarkdownドキュメントです。
H1s_HF_LDA.py の解析ドキュメント
概要
このモジュールは、H-like 1s軌道のエネルギーレベルを計算するために、1s動径関数とスレーターのX-αポテンシャルを使用します。ハートリー・フォック・スレーター（HFS）法の基本的な概念に基づき、電子間の相互作用を交換ポテンシャルで近似して扱います。動径関数の指数係数 ka や電子数 Ne などを変分法で最適化し、全エネルギーを最小化する計算も可能です。
インポートされているライブラリ
このプログラムは以下のライブラリをインポートしています。
標準ライブラリ
os: オペレーティングシステム機能へのインターフェースを提供します。sys: システム固有のパラメータと関数を提供します。math: 数学関数を提供します (exp,sqrt,log)。
非標準ライブラリ
numpy: 数値計算のための強力なN次元配列オブジェクトとツールを提供します。scipy: 科学技術計算のためのアルゴリズムとツールを提供します。scipy.integrate.quad: 数値積分機能。scipy.interpolate.interp1d: 1次元データ補間機能。scipy.optimize.minimize: 関数最小化機能。
matplotlib.pyplot: グラフ描画機能を提供します。
グローバル定数と変数
プログラムの動作を制御するための多くのグローバル定数と変数が定義されています。
物理定数
pi(3.14159265358979323846): 円周率。h(6.6260755e-34): プランク定数 (Js)。hbar(1.05459e-34): ディラック定数 (Js)。c(2.99792458e8): 光速 (m/s)。e(1.60218e-19): 素電荷 (C)。kB(1.380658e-23): ボルツマン定数 (J K^-1)。me(9.1093897e-31): 電子質量 (kg)。e0(8.854418782e-12): 真空の誘電率 (C^2 N^-1 m^-2)。e2_4pie0(2.30711e-28):e^2 / (4 * pi * e0)の値 (N m^2)。a0(5.29177e-11 * 1.0e10): ボーア半径をオングストローム単位で表した値。HartreeToeV(27.2116): ハートリー単位を電子ボルトに変換するための係数。RyToeV(HartreeToeV / 2.0): リドベリ単位を電子ボルトに変換するための係数。pi2(pi + pi):2 * pi。pi4(pi2 + pi2):4 * pi。
計算モードと出力設定
mode(初期値:'k'): 実行モードを制御する文字列。以下の文字の組み合わせで指定できます。'd': デバッグモード。基本的なグラフを表示します。'g': グラフをプロットします。'k':kaを掃引します。'n':Neを掃引します。'v': 変分計算を追加します。'e': 全エネルギー (Etot) を基準に出力します (デフォルトはE 1s固有値)。
ELabel(初期値:'E 1s'): エネルギーラベル。modeに'e'が含まれる場合は'Etot'に設定されます。
核および軌道パラメータ
Z(初期値: 1.0): 原子番号。n(初期値: 1): 主量子数 (1s軌道の場合)。l(初期値: 0): 方位量子数 (1s軌道の場合)。m(初期値: 0): 磁気量子数 (1s軌道の場合)。ka(初期値: 1.0): 1s動径関数の指数部における係数。Ne(初期値: 1.0): 電子数。alpha(初期値: 2.0 / 3.0): スレーターのX-αポテンシャルにおける交換パラメータ。
動径範囲と積分パラメータ
Rmin(初期値: 0.0): 積分範囲の最小半径。Rmax(初期値: 20.0): 積分範囲の最大半径。nR(初期値: 2001): 動径グリッドの点数。Rstep(初期値:None): 動径グリッドのステップサイズ。main()関数で計算されます。Rmaxdata(初期値:None): 動径グリッドの最大値。main()関数で計算されます。nmaxdiv(初期値: 40):scipy.integrate.quad関数の最大サブインターバル数。epsR(初期値: 1.0e-4): 動径積分における許容誤差 (具体的な用途はコードから確認できません)。eps(初期値: 1.0e-8):scipy.integrate.quad関数の相対許容誤差。
電子数 (Ne) と ka の掃引配列
hparab(初期値: 0.01):Neメッシュにおける全エネルギーの1次・2次導関数を計算するための微小変化量。Nearray:Neを掃引する際に使用される電子数の配列。kaarray:kaを掃引する際に使用されるkaの配列。
最適化パラメータ
method(初期値:"nelder-mead"):scipy.optimize.minimizeで使用される最適化アルゴリズム。選択可能なオプションは以下の通りです。"nelder-mead"(Downhill simplex)"powell"(Modified Powell)"cg"(conjugate gradient, Polak-Ribiere method)"bfgs"(BFGS法)"newton-cg"(Newton-CG)"trust-ncg"(信頼領域 Newton-CG 法)"dogleg"(信頼領域 dog-leg 法)"L-BFGS-B""TNC""COBYLA""SLSQP""trust-constr""trust-exact""trust-krylov"
maxiter(初期値: 1000): 最適化アルゴリズムの最大反復回数。tol(初期値: 1.0e-3): 最適化の許容誤差。h_diff(初期値: 1.0e-3): 1次導関数を数値的に計算する際の中心差分法における微小変化量。
動径分布データ
yRr: 動径分布関数Rr(r)の値を格納するリスト。yQr: 核からの距離rの内側に存在する総電荷Q(r)の値を格納するリスト。
関数
pfloat(str)
概要: 文字列を浮動小数点数に安全に変換します。
詳細説明: 標準の
float()関数と異なり、変換に失敗した場合でもエラーを発生させずにNoneを返します。これにより、不正な入力があってもプログラムが中断されるのを防ぎます。:param
str: 変換する文字列。:type
str: str:returns: 変換された浮動小数点数、または変換に失敗した場合は
None。:rtype: float or None
pint(str)
概要: 文字列を整数に安全に変換します。
詳細説明: 標準の
int()関数と異なり、変換に失敗した場合でもエラーを発生させずにNoneを返します。これにより、不正な入力があってもプログラムが中断されるのを防ぎます。:param
str: 変換する文字列。:type
str: str:returns: 変換された整数、または変換に失敗した場合は
None。:rtype: int or None
getarg(position, defval = None)
概要: コマンドライン引数を安全に取得します。
詳細説明:
sys.argvリストから指定された位置の引数を取得します。指定された位置が存在しない場合、エラーを発生させずにデフォルト値を返します。:param
position: 取得する引数のsys.argv内でのインデックス。:type
position: int:param
defval: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。:type
defval: any, optional:returns: 指定された位置の引数、または引数が存在しない場合は
defval。:rtype: str or any
getfloatarg(position, defval = None)
概要: コマンドライン引数を浮動小数点数として安全に取得します。
詳細説明: 指定された位置のコマンドライン引数を取得し、
pfloat()関数を使用して浮動小数点数に変換します。引数が存在しない場合や変換に失敗した場合、エラーを発生させずにデフォルト値を返します。:param
position: 取得する引数のsys.argv内でのインデックス。:type
position: int:param
defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値。デフォルトはNone。:type
defval: float or None, optional:returns: 変換された浮動小数点数、または失敗した場合は
defval。:rtype: float or None
getintarg(position, defval = None)
概要: コマンドライン引数を整数として安全に取得します。
詳細説明: 指定された位置のコマンドライン引数を取得し、
pint()関数を使用して整数に変換します。引数が存在しない場合や変換に失敗した場合、エラーを発生させずにデフォルト値を返します。:param
position: 取得する引数のsys.argv内でのインデックス。:type
position: int:param
defval: 引数が存在しない場合や変換に失敗した場合に返すデフォルト値。デフォルトはNone。:type
defval: int or None, optional:returns: 変換された整数、または失敗した場合は
defval。:rtype: int or None
usage(ka = ka, Z = Z, n = n, l = l, m = m)
概要: プログラムの使用方法を表示します。
詳細説明: コマンドライン引数の構文と利用可能なモードオプションをユーザーに示します。デフォルト値の
ka,Z,n,l,mは現在のグローバル変数から取得されますが、直接引数として渡すことも可能です。:param
ka: 1s動径関数の指数部の係数。デフォルトはグローバル変数ka。:type
ka: float, optional:param
Z: 原子番号。デフォルトはグローバル変数Z。:type
Z: float, optional:param
n: 主量子数。デフォルトはグローバル変数n。:type
n: int, optional:param
l: 方位量子数。デフォルトはグローバル変数l。:type
l: int, optional:param
m: 磁気量子数。デフォルトはグローバル変数m。:type
m: int, optional:returns: なし
:rtype: None
terminate(message = None)
概要: メッセージを表示し、プログラムを終了します。
詳細説明: オプションで終了メッセージを表示し、
usage()関数を呼び出して使用方法を提示した後、プログラムを終了します。:param
message: 終了時に表示するメッセージ。デフォルトはNone。:type
message: str, optional:returns: なし (プログラムが終了するため)
:rtype: None
updatevars()
概要: コマンドライン引数に基づいてグローバル変数を更新します。
詳細説明:
sys.argvからプログラム起動時の引数を取得し、mode,ELabel,ka,Z,Neなどのグローバル変数を更新します。これにより、ユーザーはコマンドラインから計算パラメータを指定できます。ELabelはmodeに'e'が含まれるかどうかに応じて設定されます。:returns: なし
:rtype: None
Rr(ka, Z, n, l, m, r)
概要: H-like 1s軌道の動径関数
R_nl(r)を計算します。詳細説明: この関数は、水素様原子の1s軌道に対する動径関数を計算します。正規化因子
R1s0を使用し、指数係数kaと原子番号Zに依存します。正規化条件はinteg(4pi*r*r*Rr(r)*Rr(r))[r=0, inf] = 1です。なお、この関数内で使用されるグローバル変数R1s0は、コード中では明示的に定義されていません。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数 (1s軌道なので常に1)。:type
n: int:param
l: 方位量子数 (1s軌道なので常に0)。:type
l: int:param
m: 磁気量子数 (1s軌道なので常に0)。:type
m: int:param
r: 核からの距離。:type
r: float:returns: 指定された距離におけるH-like 1s軌道の動径関数の値。
:rtype: float
rho(ka, Z, n, l, m, r)
概要: H-like 1s軌道の電子密度関数
ρ(r)を計算します。詳細説明: 動径関数
Rr(r)の二乗として電子密度を計算します。ρ(r) = |Rr(r)|^2です。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数 (1s軌道なので常に1)。:type
n: int:param
l: 方位量子数 (1s軌道なので常に0)。:type
l: int:param
m: 磁気量子数 (1s軌道なので常に0)。:type
m: int:param
r: 核からの距離。:type
r: float:returns: 指定された距離における電子密度
ρ(r)の値。:rtype: float
calQ(R)
概要: 核からの距離
Rの内側に存在する総電荷 (Q(R)) を計算します。詳細説明: 0から
Rまでの範囲で電子密度ρ(r)を積分することにより、距離Rの内側の電子が形成する電荷を求めます。この関数はbuild_Qr()で構築された補間関数qfuncを使用します。qfuncが未構築の場合や、指定されたRが有効範囲外の場合はエラー処理または0を返します。:param
R: 核からの距離。:type
R: float:returns: 距離
Rの内側に存在する総電荷。:rtype: float
build_Qr(ka, Z, n, l, m)
概要: 距離
r内側の総電荷分布Q(r)の補間関数を構築します。詳細説明: グローバル変数
rの各点に対してRr(r)とrho(r)を計算し、4 * pi * r*r * rho(r)を0からrまで積分することでQ(r)を求めます。求めたQ(r)のデータポイントyQrを使用して、scipy.interpolate.interp1dにより三次スプライン補間関数qfuncを構築し、calQ()関数で利用できるようにします。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数 (1s軌道なので常に1)。:type
n: int:param
l: 方位量子数 (1s軌道なので常に0)。:type
l: int:param
m: 磁気量子数 (1s軌道なので常に0)。:type
m: int:returns: なし。グローバル変数
yRr,yQr,qfuncを更新します。:rtype: None
calTanal(Z = 1.0, ka = 1.0)
概要: H様原子の動径関数における運動エネルギーの解析解を計算します。
詳細説明: この関数は、原子番号
Zと指数係数kaを持つH様原子の1s軌道に対する運動エネルギーの解析的な近似値をハートリー単位で計算します。T = (ka*Z)^2 / 2という関係に基づいています。:param
Z: 原子番号。デフォルトは1.0。:type
Z: float, optional:param
ka: 動径関数の指数部分の係数。デフォルトは1.0。:type
ka: float, optional:returns: 運動エネルギーの解析解 (ハートリー単位)。
:rtype: float
calUanal(Z = 1.0, ka = 1.0)
概要: H様原子の動径関数における核電子相互作用ポテンシャルエネルギーの解析解を計算します。
詳細説明: この関数は、原子番号
Zと指数係数kaを持つH様原子の1s軌道に対する核と電子の間のポテンシャルエネルギーの解析的な近似値をハートリー単位で計算します。U = -(ka*Z)^2という関係に基づいています。:param
Z: 原子番号。デフォルトは1.0。:type
Z: float, optional:param
ka: 動径関数の指数部分の係数。デフォルトは1.0。:type
ka: float, optional:returns: 核電子相互作用ポテンシャルエネルギーの解析解 (ハートリー単位)。
:rtype: float
integrate_trapezoidal(func, E0, E1, h)
概要: 台形公式を用いて数値積分を行います。
詳細説明: 指定された関数
funcをE0からE1までの範囲で、ステップサイズhを用いて台形公式で数値積分します。この関数は現在使用されていませんが、将来的な高速化のために残されています。:param
func: 積分する関数。引数を一つ取る関数を想定。:type
func: callable:param
E0: 積分範囲の下限。:type
E0: float:param
E1: 積分範囲の上限。:type
E1: float:param
h: 積分ステップサイズ。:type
h: float:returns: [積分の結果, エラー値 (ここでは常に-1.0)]
:rtype: list of float
integrate3DR(func, rmin, rmax, limit = 15, epsrel = 1.0e-8)
概要: 3次元空間における動径方向の積分 (
4 * pi * r^2 * func(r) dr) を行います。詳細説明:
scipy.integrate.quadを使用して、球座標における体積要素4 * pi * r*rを考慮した動径方向の積分を実行します。:param
func: 積分する動径関数。引数を一つ取る関数を想定。:type
func: callable:param
rmin: 積分範囲の下限。:type
rmin: float:param
rmax: 積分範囲の上限。:type
rmax: float:param
limit:quad関数の最大サブインターバル数。デフォルトは15。:type
limit: int, optional:param
epsrel:quad関数の相対許容誤差。デフォルトは1.0e-8。:type
epsrel: float, optional:returns: [積分の結果, 積分の絶対誤差]
:rtype: tuple of (float, float)
calUZ(r, Z)
概要: 核電荷
Zによって生成される静電ポテンシャルU_Z(r)を計算します。詳細説明: 点電荷としての核が距離
rに生成するクーロンポテンシャル-Z/rを計算します。rが非常に小さい場合 (1.0e-3より小さい場合) は、数値的な安定性のために-Z/1.0e-3を返します。:param
r: 核からの距離。:type
r: float:param
Z: 原子番号 (核電荷)。:type
Z: float:returns: 核電荷
Zによって生成される静電ポテンシャル。:rtype: float
calUrho(r, ka, Z, n, l, m)
概要: 電子密度
ρ(r)によって距離rに形成される静電ポテンシャルU_ee(r)を計算します。詳細説明: 電子密度
ρ(r)が自身によって生成するポテンシャルを計算します。これは以下の2つの部分から構成されます。距離
rの内側にある電荷Q(r)によるポテンシャル (Q(r)/r)。距離
rの外側にある電荷が距離rに生成するポテンシャル (4 * pi * ∫rm * ρ(rm) drm,rからinf)。
:param
r: 核からの距離。:type
r: float:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:returns: 電子密度
ρ(r)によって距離rに形成される静電ポテンシャル。:rtype: float
calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: 電子系の運動エネルギーを計算します。
詳細説明: シュレーディンガー方程式の運動エネルギー演算子
(-1/2)∇^2に基づいて、H-like 1s軌道の電子の運動エネルギーを数値積分によって計算します。積分の範囲はRminからRmaxint(実質的にはinf) です。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数 (この関数では直接使用されないが引数として存在)。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [運動エネルギー, 積分の絶対誤差]
:rtype: tuple of (float, float)
calUeZ(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: 電子と原子核間の引力ポテンシャルエネルギー
U_eZを計算します。詳細説明: 電子密度
ρ(r)と核電荷Zの間のクーロン引力によるポテンシャルエネルギーを数値積分によって計算します。積分はRminからRmaxint(実質的にはinf) です。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数 (この関数では直接使用されないが引数として存在)。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [電子と原子核間の引力ポテンシャルエネルギー, 積分の絶対誤差]
:rtype: tuple of (float, float)
calUee(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: 電子間のクーロン反発ポテンシャルエネルギー
U_eeを計算します。詳細説明: 電子密度
ρ(r)によって生成される静電ポテンシャルU_rho(r)とρ(r)の積を全空間で積分することにより、電子間のクーロン反発エネルギーを計算します。modeに'e'が含まれる場合、全エネルギー計算モードとしてNe*Ne倍、それ以外の場合はNe倍します。:param
mode: 現在の実行モードを示す文字列（'e'が含まれるかでスケールが変わる）。:type
mode: str:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [電子間のクーロン反発ポテンシャルエネルギー, 積分の絶対誤差]
:rtype: tuple of (float, float)
calLocalUXa(r, ka, Z, n, l, m, Ne)
概要: スレーターのX-α交換ポテンシャル (局所形式) を計算します。
詳細説明: 距離
rにおけるスレーターのX-α交換ポテンシャルの局所形式-3.0 * alpha * (3/(4pi) * ρ(r))^(1/3)を計算します。:param
r: 核からの距離。:type
r: float:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数 (この関数では直接使用されないが引数として存在)。:type
Ne: float:returns: 距離
rにおけるX-α交換ポテンシャルの値。:rtype: float
calUXa(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: X-α交換エネルギーを計算します。
詳細説明: スレーターのX-α交換ポテンシャル
calLocalUXa(r)と電子密度rho(r)の積を全空間で積分することにより、X-α交換エネルギーを計算します。modeに'e'が含まれる場合、全エネルギー計算モードとしてNe^(4/3)倍、それ以外の場合はNe^(1/3)倍します。:param
mode: 現在の実行モードを示す文字列（'e'が含まれるかでスケールが変わる）。:type
mode: str:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [X-α交換エネルギー, 積分の絶対誤差]
:rtype: tuple of (float, float)
calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: H-like 1s軌道の全エネルギーを計算します。
詳細説明: 運動エネルギー (
T)、核電子引力エネルギー (UeZ)、電子間反発エネルギー (Uee)、およびX-α交換エネルギー (UXa) の各項を計算し、それらを合計して全エネルギーを求めます。build_Qr()を呼び出してQ(r)分布を事前に構築します。:param
ka: 動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [運動エネルギー
T, 核電子引力エネルギーUeZ, 電子間反発エネルギーUee, X-α交換エネルギーUXa, 全エネルギーEtot]:rtype: tuple of (float, float, float, float, float)
calTotalEnergyOnly(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: 変分計算のために、特定の
ka値におけるH-like 1s軌道の全エネルギーを計算します。詳細説明: この関数は
calTotalEnergy()と同様に全エネルギーを計算しますが、最適化ルーチン (minimize) から呼び出されることを想定し、引数として与えられたkavを使用し、最終的な全エネルギー値のみを返します。また、電子間反発エネルギー (Uee) とX-α交換エネルギー (UXa) の計算モードは常に全エネルギー計算 ('e'モード) として扱われます。:param
kav: 変分パラメータとして使用される動径関数の指数部分の係数。:type
kav: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: 計算された全エネルギー。
:rtype: float
calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
概要: 変分法を用いて全エネルギーを最小化する
ka値を探索し、そのときの全エネルギーと各エネルギー成分を計算します。詳細説明:
scipy.optimize.minimize関数を使用して、動径関数の指数係数kaを変分パラメータとして、calTotalEnergyOnly()で計算される全エネルギーを最小化します。指定された最適化手法 (method) とパラメータ (tol,maxiter) に基づいて探索を行い、最小化されたkaとそのときの各エネルギー成分、および全エネルギーを返します。:param
ka: 最適化の初期値となる動径関数の指数部分の係数。:type
ka: float:param
Z: 原子番号。:type
Z: float:param
n: 主量子数。:type
n: int:param
l: 方位量子数。:type
l: int:param
m: 磁気量子数。:type
m: int:param
Ne: 電子数。:type
Ne: float:param
Rmin: 積分範囲の最小半径。:type
Rmin: float:param
Rmax: 積分範囲の最大半径（上限）。:type
Rmax: float:param
nmaxdiv:scipy.integrate.quadの最大サブインターバル数。:type
nmaxdiv: int:param
eps:scipy.integrate.quadの相対許容誤差。:type
eps: float:returns: [運動エネルギー
T, 核電子引力エネルギーUeZ, 電子間反発エネルギーUee, X-α交換エネルギーUXa, 最小化された全エネルギーEtot, 最適化されたka値]:rtype: tuple of (float, float, float, float, float, float)
diff1(ai)
概要: 1変数関数の勾配 (1次導関数) を中心差分法で近似します。
詳細説明:
scipy.optimize.minimizeの勾配情報が必要なメソッド (cgなど) のために、パラメータai(ここではkaの配列) の微小変化h_diffを用いて、calTotalEnergyOnly()関数に対する1次導関数を数値的に計算します。この実装では、入力aiが変更される可能性があります。:param
ai: 最適化される変数の現在の値 (ここではkaの配列)。:type
ai: numpy.ndarray:returns: 各変数に対する1次導関数の配列。
:rtype: numpy.ndarray
callback(xk)
概要: 最適化の各イテレーションで呼び出され、進行状況を表示します。
詳細説明:
scipy.optimize.minimize関数のcallbackオプションとして使用されます。各イテレーションで現在のパラメータxk(ここではkaの配列) とそれに対応する全エネルギーEtotを表示します。:param
xk: 現在のイテレーションでの最適化変数の値 (ここではkaの配列)。:type
xk: numpy.ndarray:returns: なし
:rtype: None
sweep_Ne(Eanal)
概要: 電子数
Neを掃引し、非最適化および最適化された全エネルギーの変化を評価・プロットします。詳細説明:
Nearrayに定義された様々な電子数Neに対して、固定されたkaでの全エネルギー、および変分法でkaを最適化した上での全エネルギーを計算します。また、Ne=0.5付近での放物線近似も行い、結果を表形式で出力し、グラフとしてプロットします。特に 'v' モードが有効な場合、kaの最適化も行います。グラフの表示は 'g' モードが有効な場合のみ行われます。:param
Eanal: 比較のための解析解エネルギー値。:type
Eanal: float:returns: なし
:rtype: None
sweep_ka(Eanal)
概要: 動径関数の指数係数
kaを掃引し、全エネルギーの変化を評価・プロットします。詳細説明:
kaarrayに定義された様々なkaの値に対して、H-like 1s軌道の全エネルギーとその各成分を計算します。結果を表形式で出力し、'g' モードが有効な場合はグラフとしてプロットします。:param
Eanal: 比較のための解析解エネルギー値。:type
Eanal: float:returns: なし
:rtype: None
debug(Eanal)
概要: デバッグモードで、H-like 1s軌道の計算結果と関連するグラフを表示します。
詳細説明: 現在のパラメータ (
ka,Z,Neなど) における運動エネルギー、核電子引力エネルギー、電子間反発エネルギー、X-α交換エネルギー、および全エネルギーを計算し、表形式で出力します。'v' モードが有効な場合は、kaを最適化した結果も表示します。また、'g' モードが有効な場合は、動径関数Rr(r)、電荷分布Q(r)、および各種ポテンシャル (U(Z),U(Z-rho),U(rho),U(Xa)) のグラフをプロットします。:param
Eanal: 比較のための解析解エネルギー値。:type
Eanal: float:returns: なし
:rtype: None
main()
概要: プログラムのメインエントリポイントです。
詳細説明: コマンドライン引数を解析し、グローバル変数を更新します。水素様原子の解析解エネルギーを計算した後、指定されたモードに基づいて異なる計算（デバッグ、
ka掃引、Ne掃引）を実行します。最後に、動径関数の正規化チェックも行われます。グローバル変数rとRmaxdataはこの関数で初期化されます。:returns: なし
:rtype: None
使用方法
H1s_HF_LDA.py はコマンドライン引数を受け取ります。
python H1s_HF_LDA.py mode Z ka Ne
mode: 以下のキー文字を組み合わせた文字列。d: デバッグモード。基本的なグラフをプロットします。v: 変分計算を追加します。e: 全エネルギー (Etot) を基準に出力します (デフォルトはE 1s固有値)。k:kaを掃引します。n:Neを掃引します。g: グラフをプロットします。
Z: 原子番号 (浮動小数点数)。ka: 1s動径関数の指数部の係数 (浮動小数点数)。Ne: 電子数 (浮動小数点数)。
引数はオプションであり、省略された場合はグローバル変数のデフォルト値が使用されます。
例:
python H1s_HF_LDA.py nvg 1.0 1.0 0.5
python H1s_HF_LDA.py k 1.0 1.0 1.0
入出力
入力
プログラムはコマンドライン引数から以下のパラメータを受け取ります。
計算モード (
mode): 実行する計算の種類 (例:'d','k','n','v','g','e'の組み合わせ)。原子番号 (
Z): 計算対象の原子の核電荷。動径関数の指数係数 (
ka): 1s軌道動径関数exp(-ka * Z * r)のka値。電子数 (
Ne): 計算に含める電子数。
これらの引数は updatevars() 関数によって解析され、対応するグローバル変数が設定されます。
出力
標準出力
プログラム設定: 起動時に
mode, 軌道パラメータ (ka,Z,n,l,m), 電子数 (Ne), 積分範囲 (Rmax) など、現在の計算設定が表示されます。解析解: 水素様原子の運動エネルギー (
Tanal)、ポテンシャルエネルギー (Uanal)、および全エネルギー (Eanal) の解析解が電子ボルト (eV) 単位で表示されます。正規化チェック: 動径関数の正規化積分の結果 (
2 * pi * integ(r*r * Rr*2)dr) が表示され、値が1に近似するかどうかを確認できます。計算結果の表:
現在の
ka,Z,Neにおける、運動エネルギー (T), 核電子引力エネルギー (UeZ), 電子間反発エネルギー (Uee), X-α交換エネルギー (UXa), および全エネルギー (EtotまたはE 1s) が電子ボルト (eV) 単位で表形式で出力されます。変分法による最適化が有効な場合 (
'v'モード)、最適化されたka値とそれに対応するエネルギー成分も出力されます。sweep_Ne()関数が実行される場合 ('n'モード)、Ne=0.5周辺でのエネルギーの放物線近似の結果と、Nearrayの各Ne値に対する計算結果が表形式で出力されます。sweep_ka()関数が実行される場合 ('k'モード)、kaarrayの各ka値に対する計算結果が表形式で出力されます。
最適化の進捗:
scipy.optimize.minimizeのコールバック関数 (callback()) により、最適化の各イテレーションにおける現在のka値と全エネルギーがコンソールに表示されます。終了プロンプト: グラフが表示された場合、ユーザーがEnterキーを押すまでプログラムは待機します。
グラフィック出力
'g' モードが指定されている場合、matplotlib.pyplot を使用して以下のグラフが表示されます。
debug()モード ('d'および'g'):動径関数: 核からの距離
rに対する動径関数Rr(r)のグラフ。総電荷分布: 核からの距離
rの内側に存在する総電荷Q(r)のグラフ。各種ポテンシャル: 核からの距離
rに対する核ポテンシャル (U(Z))、電子-核ポテンシャル (U(Z-rho))、電子間ポテンシャル (-U(rho))、およびX-α交換ポテンシャル (U(Xa)) のグラフ。
sweep_Ne()モード ('n'および'g'):エネルギー対
Ne: 電子数Neに対する全エネルギー (ELabel) の変化のグラフ。非最適化の結果と、'v' モードが有効な場合は最適化された結果、およびそれぞれの放物線近似が表示されます。最適化された
ka対Ne: 電子数Neに対する最適化されたka値の変化のグラフ ('v'モードが有効な場合のみ)。
sweep_ka()モード ('k'および'g'):エネルギー対
ka: 動径関数の指数係数kaに対する全エネルギー (ELabel) の変化のグラフ。