"""
H-like 1s軌道エネルギーレベル計算モジュール
概要:
このモジュールは、H-like 1s軌道のエネルギーレベルを計算するために、1s動径関数とスレーターのX-αポテンシャルを使用します。
ハートリー・フォック・スレーター（HFS）法の基本的な概念に基づき、電子間の相互作用を交換ポテンシャルで近似して扱います。
動径関数の指数係数kaや電子数Neなどを変分法で最適化し、全エネルギーを最小化する計算も可能です。
関連リンク:
:doc:`H1s-HF-LDA_usage`
"""
import os
import sys
from math import exp, sqrt
import numpy as np
from math import log, exp
from numpy import arange
from scipy import integrate # 数値積分関数 integrateを読み込む
from scipy.interpolate import interp1d
from scipy import optimize
from scipy.optimize import minimize
from matplotlib import pyplot as plt
# constants
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 # JK<sup>-1</sup>";
me = 9.1093897e-31 # kg";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e2_4pie0 = 2.30711e-28 # Nm<sup>2</sup>";
a0 = 5.29177e-11 * 1.0e10 # A
HartreeToeV = 27.2116 # eV";
RyToeV = HartreeToeV / 2.0
pi2 = pi + pi
pi4 = pi2 + pi2
# mode: 'd': debug mode, plot fundamental graphs
# 'g': plot graph
# 'k': sweep ka, 'n': sweep Ne
# 'v': add variational calculation, 'e': energy-based output
mode = 'k'
ELabel = 'E 1s'
# Nuclear and orbital parameters
Z = 1.0
n = 1
l = 0
m = 0
ka = 1.0 # coefficient of the exponent in R1s
Ne = 1.0
# Number of electrons, Slater's alpha
alpha = 2.0 / 3.0
# Radius range and integration (quad()) parameters
Rmin = 0.0
Rmax = 20.0
nR = 2001
Rstep = None
Rmaxdata = None
nmaxdiv = 40
epsR = 1.0e-4
eps = 1.0e-8
# Ne mesh to calculate 1st and 2nd derivative of Etot
hparab = 0.01
# iteration parameters
#Nearray = [1.0, 0.8, 0.6]
#Nearray = [1.0, 0.8, 0.6, 0.5, 0.4, 0.3, 0.1, 0.01, 0.001]
Nearray = [1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.07, 0.05, 0.04, 0.03, 0.02, 0.01, 1.0e-3, 1.0e-4, 1.0e-5, 1.0e-6]
kaarray = [0.5, 0.6, 0.65, 0.7, 0.8, 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4]
#=============================================
# 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’ (see here)
#TNC’ (see here)
#COBYLA’ (see here)
#SLSQP’ (see here)
#trust-constr’(see here)
#dogleg’ (see here)
#trust-exact’ (see here)
#trust-krylov’ (see here)
method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'
maxiter = 1000
tol = 1.0e-3
h_diff = 1.0e-3
[ドキュメント]
def pfloat(str):
"""
概要: 文字列を浮動小数点数に安全に変換します。
詳細説明:
標準のfloat()関数と異なり、変換に失敗した場合でもエラーを発生させずにNoneを返します。
これにより、不正な入力があってもプログラムが中断されるのを防ぎます。
:param str: 変換する文字列。
:type str: str
:returns: 変換された浮動小数点数、または変換に失敗した場合はNone。
:rtype: float or None
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
概要: 文字列を整数に安全に変換します。
詳細説明:
標準のint()関数と異なり、変換に失敗した場合でもエラーを発生させずにNoneを返します。
これにより、不正な入力があってもプログラムが中断されるのを防ぎます。
:param str: 変換する文字列。
:type str: str
:returns: 変換された整数、または変換に失敗した場合はNone。
:rtype: int or None
"""
try:
return int(str)
except:
return None
[ドキュメント]
def 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
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def 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
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def 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
"""
return pint(getarg(position, defval))
[ドキュメント]
def 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
"""
print("")
print("Usage: Variables in () are optional")
print(" (i) python {} mode Z ka Ne".format(sys.argv[0]))
print(" mode: combination of the following key characters")
print(" d: debug mode, plot fundamental graphs")
print(" v: add variational calculations")
print(" e: output based on energy (default: based on E 1s eigen value)")
print(" k: sweep ka")
print(" n: sweep Ne")
print(" g : plot graph")
print(" ex: python {} nvg {} {} {}".format(sys.argv[0], Z, ka, Ne))
print(" ex: python {} k {} {} {}".format(sys.argv[0], Z, ka, Ne))
[ドキュメント]
def terminate(message = None):
"""
概要: メッセージを表示し、プログラムを終了します。
詳細説明:
オプションで終了メッセージを表示し、usage()関数を呼び出して使用方法を提示した後、
プログラムを終了します。
:param message: 終了時に表示するメッセージ。デフォルトはNone。
:type message: str, optional
:returns: なし (プログラムが終了するため)
:rtype: None
"""
if message is not None:
print("")
print(message)
usage()
print("")
exit()
[ドキュメント]
def updatevars():
"""
概要: コマンドライン引数に基づいてグローバル変数を更新します。
詳細説明:
sys.argvからプログラム起動時の引数を取得し、`mode`, `ELabel`, `ka`, `Z`, `Ne`
などのグローバル変数を更新します。これにより、ユーザーはコマンドラインから計算パラメータを
指定できます。`ELabel`は`mode`に'e'が含まれるかどうかに応じて設定されます。
:returns: なし
:rtype: None
"""
global mode, ELabel
global ka, Z, Ne, n, l, m, Ne, Rmax, Rstep
argv = sys.argv
if len(argv) == 1:
terminate()
mode = getarg( 1, mode)
Z = getfloatarg( 2, Z)
ka = getfloatarg( 3, ka)
Ne = getfloatarg( 4, Ne)
if 'e' in mode:
ELabel = 'Etot'
else:
ELabel = 'E 1s'
# Total charge inside r
yRr = [] # Raidal distributuion function
yQr = [] # Integrated charge inside r: integ_rm(4pi * rm * rm * rho(rm))_[rm = 0, r]
[ドキュメント]
def 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 です。
: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
"""
global R1s0
return R1s0 * pow(ka * Z, 1.5) * exp(-ka * Z * r)
[ドキュメント]
def 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
"""
psi = Rr(ka, Z, n, l, m, r)
return psi * psi
[ドキュメント]
def calQ(R):
"""
概要: 核からの距離Rの内側に存在する総電荷 (Q(R)) を計算します。
詳細説明:
0からRまでの範囲で電子密度 ρ(r) を積分することにより、
距離Rの内側の電子が形成する電荷を求めます。
この関数は `build_Qr` で構築された補間関数 `qfunc` を使用します。
`qfunc`が未構築の場合や、指定されたRが有効範囲外の場合はエラー処理または0を返します。
:param R: 核からの距離。
:type R: float
:returns: 距離Rの内側に存在する総電荷。
:rtype: float
"""
global Rmin, Rmaxdata, r, Qr
if R < Rmin:
print("Error in Q(r): Given r={} exceeds the R range [{}, {}]".format(R, Rmin, Rmax))
exit()
if R > Rmaxdata:
return 0.0
return qfunc(R)
[ドキュメント]
def build_Qr(ka, Z, n, l, m):
"""
概要: 距離r内側の総電荷分布 Q(r) の補間関数を構築します。
詳細説明:
`r`の各点に対して `Rr(r)` と `rho(r)` を計算し、
`4pi * 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
"""
global yRr, yQr, qfunc
yRr = []
yQr = []
for i in range(len(r)):
yRr.append(Rr(ka, Z, n, l, m, r[i]))
if r[i] == 0.0:
Q = 0.0
else:
Q, errQ = integrate.quad(lambda r: pi4 * r * r * rho(ka, Z, n, l, m, r), 0, r[i], limit = nmaxdiv, epsrel = eps)
yQr.append(Ne * Q)
qfunc = interp1d(r, yQr, kind = 'cubic')
[ドキュメント]
def 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
"""
# T = e^2 / 8pi / e0 / a0
Tanal = Z * Z * ka * ka * e / 8.0 / pi / e0 / (a0*1.0e-10)
return Tanal
[ドキュメント]
def 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
"""
# U = -e^2 / 4pi / e0 / a0
Uanal = -Z * Z * ka * ka * e / 4.0 / pi / e0 / (a0*1.0e-10)
return Uanal
[ドキュメント]
def 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
"""
n = int((E1 - E0) / h + 1.000001)
h = (E1 - E0) / n
y = [func(E0 + i * h) for i in range(n)]
S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1])
return [h * S, -1.0]
[ドキュメント]
def integrate3DR(func, rmin, rmax, limit = 15, epsrel = 1.0e-8):
"""
概要: 3次元空間における動径方向の積分 (4πr^2 func(r) dr) を行います。
詳細説明:
`scipy.integrate.quad` を使用して、球座標における体積要素 `4pi * 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)
"""
I, err = integrate.quad(lambda r: pi4 * r * r * func(r), rmin, rmax, limit, epsrel)
return I, err
[ドキュメント]
def 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
"""
if r < 1.0e-3:
return -Z / 1.0e-3
else:
return -Z / r
[ドキュメント]
def calUrho(r, ka, Z, n, l, m):
"""
概要: 電子密度 ρ(r) によって距離rに形成される静電ポテンシャル U_ee(r) を計算します。
詳細説明:
電子密度 ρ(r) が自身によって生成するポテンシャルを計算します。
これは以下の2つの部分から構成されます。
1. 距離rの内側にある電荷 Q(r) によるポテンシャル (Q(r)/r)。
2. 距離rの外側にある電荷が距離rに生成するポテンシャル (4pi * ∫rm * ρ(rm) drm, rから∞)。
: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
"""
# Uerho(r) = integ(rho(m) / |r-rm|)_[rm=0, inf] [r, rm are vectors]
# = Q(r) / r + integ[4pi * rm * rm * rho(rm) / rm]_[rm = r, inf]) [r is vector]
if r == 0.0:
Uee1 = 0.0
else:
Uee1 = calQ(r) / r
Rmaxint = Rmax
# Rmaxint = min(-log(eps) / Z / ka, Rmax)
Uee2, errUee2 = integrate.quad(lambda rm: pi4 * rm * rho(ka, Z, n, l, m, rm), r, Rmaxint, limit = nmaxdiv, epsrel = eps)
return Uee1 + Uee2
# return Uee1 + Ne * Uee2
[ドキュメント]
def calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
"""
概要: 電子系の運動エネルギーを計算します。
詳細説明:
シュレーディンガー方程式の運動エネルギー演算子 `(-1/2)∇^2` に基づいて、
H-like 1s軌道の電子の運動エネルギーを数値積分によって計算します。
積分の範囲は `Rmin` から `Rmaxint` (実質的には∞) です。
: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)
"""
Rmaxint = min(-log(eps) / Z / ka, Rmax)
T, errT = integrate.quad(lambda r: r * (2.0 - Z * ka * r) * exp(-2.0 * Z * ka * r),
Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
T *= 1.0 / 2.0 * pi4 * Z * ka * R1s0 * R1s0 * pow(ka * Z, 3.0)
return T, errT
[ドキュメント]
def calUeZ(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
"""
概要: 電子と原子核間の引力ポテンシャルエネルギー U_eZ を計算します。
詳細説明:
電子密度 ρ(r) と核電荷 Z の間のクーロン引力によるポテンシャルエネルギーを
数値積分によって計算します。積分は `Rmin` から `Rmaxint` (実質的には∞) です。
: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)
"""
Rmaxint = min(-log(eps) / Z / ka, Rmax)
UeZ, errUeZ = integrate.quad(lambda r: pi4 * r * rho(ka, Z, n, l, m, r),
Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
UeZ *= -Z
return UeZ, errUeZ
[ドキュメント]
def 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)
"""
Rmaxint = min(-log(eps) / Z / ka, Rmax)
Uee, errUee = integrate.quad(lambda r: pi4 * r * r * calUrho(r, ka, Z, n, l, m) * rho(ka, Z, n, l, m, r),
Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
if 'e' in mode:
Uee *= Ne * Ne
else:
Uee *= Ne
return Uee, errUee
[ドキュメント]
def 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
"""
return -3.0 * alpha * pow(3.0/pi4 * rho(ka, Z, n, l, m, r), 1.0/3.0)
[ドキュメント]
def 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)
"""
Rmaxint = min(-log(eps) / Z / ka, Rmax)
UXa, errUXa = integrate.quad(lambda r: pi4 * r * r * calLocalUXa(r, ka, Z, n, l, m, Ne) * rho(ka, Z, n, l, m, r),
Rmin, Rmaxint, limit = nmaxdiv, epsrel = eps)
if 'e' in mode:
UXa *= pow(Ne, 4.0/3.0)
else:
UXa *= pow(Ne, 1.0/3.0)
return UXa, errUXa
[ドキュメント]
def 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)
"""
build_Qr(ka, Z, n, l, m)
#def calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps):
T, errT = calT(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
UeZ, errUeZ = calUeZ(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Uee, errUee = calUee(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
UXa, errUXa = calUXa(mode, ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Etot = T + UeZ + Uee + UXa
return T, UeZ, Uee, UXa, Etot
[ドキュメント]
def 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
"""
build_Qr(kav, Z, n, l, m)
T, errT = calT(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
UeZ, errUeZ = calUeZ(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Uee, errUee = calUee('e', kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
UXa, errUXa = calUXa('e', kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Etot = T + UeZ + Uee + UXa
return Etot
[ドキュメント]
def 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)
"""
xa = [ka]
ret = minimize(lambda xa: calTotalEnergyOnly(xa[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps),
xa, method = method, jac = diff1, tol = tol,
callback = lambda xa: callback(xa),
options = {'maxiter':maxiter, "disp":True})
# print("ret=", ret)
if method == 'nelder-mead':
simplex = ret['final_simplex']
ai = simplex[0][0]
Emin = ret['fun']
elif method == 'cg':
ai = ret['x']
Emin = ret['fun']
elif method == 'powell':
ai = ret['x']
Emin = ret['fun']
elif method == 'bfgs':
ai = ret['x']
ka = xa[0]
Emin = calTotalEnergyOnly(xa[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
ka = ai[0]
T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
return T, UeZ, Uee, UXa, Etot, ka
[ドキュメント]
def diff1(ai):
"""
概要: 1変数関数の勾配 (1次導関数) を中心差分法で近似します。
詳細説明:
`scipy.optimize.minimize` の勾配情報が必要なメソッド (`cg`など) のために、
パラメータ `ai` (ここでは `ka`) の微小変化 `h_diff` を用いて、
`calTotalEnergyOnly` 関数に対する1次導関数を数値的に計算します。
:param ai: 最適化される変数の現在の値 (ここでは `ka` の配列)。
:type ai: numpy.ndarray
:returns: 各変数に対する1次導関数の配列。
:rtype: numpy.ndarray
"""
global Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps
n = len(ai)
f0 = calTotalEnergyOnly(ai[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
df = np.empty(n)
for i in range(n):
aii = ai
aii[i] = ai[i] + h_diff
f1 = calTotalEnergyOnly(aii[0], Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
df[i] = (f1 - f0) / h_diff
return df
iter = 0
[ドキュメント]
def callback(xk):
"""
概要: 最適化の各イテレーションで呼び出され、進行状況を表示します。
詳細説明:
`scipy.optimize.minimize` 関数の `callback` オプションとして使用されます。
各イテレーションで現在のパラメータ `xk` (ここでは `ka`) と
それに対応する全エネルギー `Etot` を表示します。
:param xk: 現在のイテレーションでの最適化変数の値 (ここでは `ka` の配列)。
:type xk: numpy.ndarray
:returns: なし
:rtype: None
"""
global iter
global Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps
ka = xk[0]
Etot = calTotalEnergyOnly(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
# パラメータと残差二乗和を表示
print("callback {}: ka={:14.4g} Emin={} eV".format(iter, ka, Etot * HartreeToeV))
iter += 1
[ドキュメント]
def sweep_Ne(Eanal):
"""
概要: 電子数Neを掃引し、非最適化および最適化された全エネルギーの変化を評価・プロットします。
詳細説明:
`Nearray` に定義された様々な電子数Neに対して、固定された `ka` での全エネルギー、
および変分法で `ka` を最適化した上での全エネルギーを計算します。
また、Ne=0.5付近での放物線近似も行い、結果を表形式で出力し、グラフとしてプロットします。
特に 'v' モードが有効な場合、`ka` の最適化も行います。
:param Eanal: 比較のための解析解エネルギー値。
:type Eanal: float
:returns: なし
:rtype: None
"""
global Rmin, Rstep, nR, Rmaxdata, r
global ka
global qfunc
print("")
print("Not optimized")
print(" Parabolic expantion around Ne = 0.5")
Tp, UeZp, Ueep, UXap, Etotp = calTotalEnergy(ka, Z, n, l, m, 0.5 + hparab, Rmin, Rmax, nmaxdiv, eps)
T0, UeZ0, Uee0, UXa0, Etot0 = calTotalEnergy(ka, Z, n, l, m, 0.5 , Rmin, Rmax, nmaxdiv, eps)
Tm, UeZm, Ueem, UXam, Etotm = calTotalEnergy(ka, Z, n, l, m, 0.5 - hparab, Rmin, Rmax, nmaxdiv, eps)
Etotp *= HartreeToeV
Etot0 *= HartreeToeV
Etotm *= HartreeToeV
print(" Ne={}: E={} eV".format(0.5 + hparab, Etotp))
print(" Ne={}: E={} eV".format(0.5 , Etot0))
print(" Ne={}: E={} eV".format(0.5 - hparab, Etotm))
a2 = (Etotp - 2.0 * Etot0 + Etotm) / hparab / hparab
a1 = (Etotp - Etotm) / 2.0 / hparab
print(" E = {} + {} * Ne + {} * Ne^2".format(Etot0, a1, a2))
xNe = []
yE1s = []
yE1sparab = []
print("")
print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)',
ELabel, ELabel + ' (parbolic)'))
for Ne in Nearray:
T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Eparab = Etot0 + a1 * (Ne - 0.5) + a2 * pow(Ne - 0.5, 2)
print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
.format(ka, Z, Ne,
T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV,
Etot * HartreeToeV, Eparab))
xNe.append(Ne)
yE1s.append(Etot * HartreeToeV)
yE1sparab.append(Eparab)
if 'v' in mode:
yE1sOpt = []
yE1sOptparab = []
ykaOpt = []
print("")
print("Optimized for ka by variational principle")
print(" Parabolic expantion around Ne = 0.5")
Tp, UeZp, Ueep, UXap, Etotp, kap = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5 + hparab, Rmin, Rmax, nmaxdiv, eps)
T0, UeZ0, Uee0, UXa0, Etot0, ka0 = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5 , Rmin, Rmax, nmaxdiv, eps)
Tm, UeZm, Ueem, UXam, Etotm, kam = calOptimizedTotalEnergy(ka, Z, n, l, m, 0.5 - hparab, Rmin, Rmax, nmaxdiv, eps)
Etotp *= HartreeToeV
Etot0 *= HartreeToeV
Etotm *= HartreeToeV
print(" Ne={}: E={} eV".format(0.5 + hparab, Etotp))
print(" Ne={}: E={} eV".format(0.5 , Etot0))
print(" Ne={}: E={} eV".format(0.5 - hparab, Etotm))
a2 = (Etotp - 2.0 * Etot0 + Etotm) / hparab / hparab
a1 = (Etotp - Etotm) / 2.0 / hparab
print(" E = {} + {} * Ne + {} * Ne^2".format(Etot0, a1, a2))
print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel, ELabel + ' (parabolic)'))
for Ne in Nearray:
T, UeZ, Uee, UXa, Etot, ka = calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
Eparab = Etot0 + a1 * (Ne - 0.5) + a2 * pow(Ne - 0.5, 2)
print("{:6.4f}\t{:6.4f}\t{:6.4g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}t{:10.6g}"
.format(ka, Z, Ne,
T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV,
Etot * HartreeToeV, Eparab))
yE1sOpt.append(Etot * HartreeToeV)
yE1sOptparab.append(Eparab)
ykaOpt.append(ka)
#=============================
# Plot graphs
#=============================
if 'g' not in mode:
terminate()
fig = plt.figure(figsize = (12, 8))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
# ax3 = fig.add_subplot(2, 3, 3)
# ax4 = fig.add_subplot(2, 3, 4)
# ax5 = fig.add_subplot(2, 3, 5)
ax1.plot(xNe, yE1s, label = ELabel + ' (non-optimized)', color = 'black', linewidth = 1.5, marker = 'o', markersize = 1.0)
ax1.plot(xNe, yE1sparab, label = ELabel + ' (non-opt,parabolic)', color = 'black', linestyle = 'dashed', linewidth = 0.5)
if 'v' in mode:
ax1.plot(xNe, yE1sOpt, label = ELabel + ' (optimized)', color = 'red', linewidth = 1.5, marker = 'o', markersize = 1.0)
ax1.plot(xNe, yE1sOptparab, label = ELabel + ' (opt,parabolic)', color = 'red', linewidth = 0.5)
ax1.plot([min(xNe), max(xNe)], [Eanal, Eanal], color = 'red', linestyle = 'dashed')
ax1.plot([0.5, 0.5], [min(yE1s), max(yE1s)], color = 'red', linestyle = 'dashed')
ax1.set_xlabel("Ne")
ax1.set_ylabel(ELabel + " (eV)")
ax1.legend()
if 'v' in mode:
ax2.plot(xNe, ykaOpt, label = 'ka (optimized)', color = 'black', linewidth = 0.5, marker = 'o', markersize = 1.0)
ax2.plot([min(xNe), max(xNe)], [1.0, 1.0], color = 'red', linestyle = 'dashed')
if 'v' in mode:
ax2.plot([0.5, 0.5], [min(ykaOpt), max(ykaOpt)], color = 'red', linestyle = 'dashed')
ax2.set_xlabel("Ne")
ax2.set_ylabel("ka (optimized)")
ax2.legend()
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def sweep_ka(Eanal):
"""
概要: 動径関数の指数係数kaを掃引し、全エネルギーの変化を評価・プロットします。
詳細説明:
`kaarray` に定義された様々な `ka` の値に対して、
H-like 1s軌道の全エネルギーとその各成分を計算します。
結果を表形式で出力し、'g' モードが有効な場合はグラフとしてプロットします。
:param Eanal: 比較のための解析解エネルギー値。
:type Eanal: float
:returns: なし
:rtype: None
"""
global Rmin, Rstep, nR, Rmaxdata, r
global qfunc
global nmaxdiv, eps
xka = []
yE1s = []
print("")
print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
for kav in kaarray:
T, UeZ, Uee, UXa, Etot = calTotalEnergy(kav, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
.format(kav, Z, Ne,
T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))
xka.append(kav)
yE1s.append(Etot * HartreeToeV)
#=============================
# Plot graphs
#=============================
if 'g' not in mode:
terminate()
fig = plt.figure(figsize = (12, 8))
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 3)
ax4 = fig.add_subplot(2, 3, 4)
ax5 = fig.add_subplot(2, 3, 5)
ax1.plot(xka, yE1s, label = ELabel, color = 'black', marker = 'o')
ax1.plot([min(xka), max(xka)], [Eanal, Eanal], color = 'red', linestyle = 'dashed')
ax1.set_xlabel("ka")
ax1.set_ylabel(ELabel + " (eV)")
ax1.legend()
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def 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
"""
global ka, Z, n, l, m, Ne
global Rmin, Rstep, nR, Rmaxdata, r
ka0, Z0 = ka, Z
print("")
print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
T, UeZ, Uee, UXa, Etot = calTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
.format(ka, Z, Ne,
T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))
if 'v' in mode:
print("")
print("Optimized for ka by variational principle")
print("{:6s}\t{:6s}\t{:6s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('ka', 'Z', 'Ne', 'T in eV', 'U(e-Z)', 'U(e-e)', 'U(Xa)', ELabel))
T, UeZ, Uee, UXa, Etot, ka = calOptimizedTotalEnergy(ka, Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps)
print("{:6.4f}\t{:6.4f}\t{:6.4f}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}\t{:10.6g}"
.format(ka, Z, Ne,
T * HartreeToeV, UeZ * HartreeToeV, Uee * HartreeToeV, UXa * HartreeToeV, Etot * HartreeToeV))
#=============================
# Plot graphs
#=============================
if 'g' not in mode:
terminate()
build_Qr(ka0, Z, n, l, m)
ypotZ = [] # U(Z) = -Z/r
ypotZrho = [] # U(e-Z) for Z = 1
ypotXa = [] # Xa potential
ypotrho = [] # Electrostatic potential by rho(r)
for i in range(len(r)):
if r[i] == 0.0:
phi = 0.0
else:
phi, errpot = integrate.quad(lambda r: pi4 * r * rho(ka, Z, n, l, m, r),
Rmin, r[i], limit = nmaxdiv, epsrel = eps)
potZ = -Z / r[i]
ypotZ.append(calUZ(r[i], Z))
ypotZrho.append(-Z * Ne * phi)
ypotrho.append(Ne * calUrho(r[i], ka, Z, n, l, m))
ypotXa.append(calLocalUXa(r[i], ka, Z, n, l, m, 1.0))
fig = plt.figure(figsize = (12, 8))
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 3)
ax1.plot(r, yRr, label = 'Rr(r)', color = 'black')
ax1.set_xlim([Rmin, Rmax])
ax1.set_ylim([0, max(yRr)*1.1])
ax1.set_xlabel("r (bohr)")
ax1.set_ylabel("Rr(r)")
ax1.legend()
ax2.plot(r, yQr, label = 'Q(r)', color = 'black')
ax2.set_xlabel("r (bohr)")
ax2.set_ylabel("Q(r)")
ax2.legend()
ax3.plot(r, ypotZ, label = 'U(Z)', color = 'black')
ax3.plot(r, ypotZrho, label = 'U(Z-rho)', color = 'red')
ax3.plot(r, -np.array(ypotrho), label = '-U(rho)', color = 'blue')
ax3.plot(r, ypotXa, label = 'U(Xa)', color = 'green')
ax3.set_ylim([min(ypotXa) * 5.0, 0.0])
ax3.set_xlabel("r (bohr)")
ax3.set_ylabel("Potential / Energy (Hartree)")
ax3.legend()
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def main():
"""
概要: プログラムのメインエントリポイントです。
詳細説明:
コマンドライン引数を解析し、グローバル変数を更新します。
水素様原子の解析解エネルギーを計算した後、
指定されたモードに基づいて異なる計算（デバッグ、ka掃引、Ne掃引）を実行します。
最後に、動径関数の正規化チェックも行います。
:returns: なし
:rtype: None
"""
global mode
global ka, Z, n, l, m, Ne
global Rmin, Rstep, nR, Rmaxdata, r
global qfunc
updatevars()
Rstep = (Rmax - Rmin) / (nR - 1)
r = [Rmin + i * Rstep for i in range(nR+100)]
Rmaxdata = max(r)
print("")
print("mode: ", mode)
print("")
print("Orbital: ka={} Z={} n={} l={} m={}".format(ka, Z, n, l, m))
print("Ne: ", Ne)
print("Integration: Rmax=", Rmax)
print(" Rmax: epsR={} Rmaxinteg={}".format(epsR, -log(epsR) / Z / ka))
print("")
print("Analytical solution")
Tanal = calTanal(Z, ka)
Uanal = calUanal(Z, ka)
Eanal = Tanal + Uanal
print("T(analytical) = {} eV".format(Tanal))
print("U(analytical) = {} eV".format(Uanal))
print(" Etotl(analytical) = {} eV".format(Eanal))
print("")
print("Numerical integration")
Rr2tot, err = integrate.quad(lambda r: r * r * rho(ka, Z, n, l, m, r), Rmin, Rmax, limit = nmaxdiv, epsrel = eps)
Rr2tot *= pi4
print("R(r) normalization check: 2pi * integ(r*r * Rr*2)dr = ", Rr2tot)
if 'd' in mode:
debug(Eanal)
if 'k' in mode:
sweep_ka(Eanal)
if 'n' in mode:
sweep_Ne(Eanal)
terminate()
if __name__ == '__main__':
usage()
main()