"""
Calculate energy level of H-like 1s orbital
using 1s radial function and Slator's X-alpha potential
概要:
H-like 1s軌道のエネルギー準位を、1s動径関数とSlaterのX-alphaポテンシャルを用いて計算します。
詳細説明:
このスクリプトは、変分法と数値積分を組み合わせて、H-like 1s軌道の運動エネルギー、核-電子間ポテンシャルエネルギー、
電子-電子間反発エネルギー、およびX-alpha交換エネルギーを計算します。
コマンドライン引数により、デバッグモード、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-1";
me = 9.1093897e-31 # kg";
e0 = 8.854418782e-12; # C2N-1m-2";
e2_4pie0 = 2.30711e-28 # Nm2";
a0 = 5.29177e-11 * 1.0e10 # A
HartreeToeV = 27.2116 # eV";
RyToeV = HartreeToeV / 2.0
pi2 = pi + pi
pi4 = pi2 + pi2
R1s0 = a0
# 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: (str) 変換する文字列。
:returns: (float or None) 変換された浮動小数点数、または変換できなかった場合はNone。
"""
try:
return float(str)
except:
return None
def pint(str):
"""
概要:
文字列を整数に安全に変換します。
詳細説明:
文字列をint型に変換しようと試みます。変換に失敗した場合でもプログラムを終了せず、Noneを返します。
:param str: (str) 変換する文字列。
:returns: (int or None) 変換された整数、または変換できなかった場合はNone。
"""
try:
return int(str)
except:
return None
def getarg(position, defval = None):
"""
概要:
コマンドライン引数を安全に取得します。
詳細説明:
sys.argvリストから指定された位置の引数を取得します。
指定された位置に引数が存在しない場合、エラーを発生させずにデフォルト値を返します。
:param position: (int) 取得する引数の位置（インデックス）。
:param defval: (any, optional) 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
:returns: (str or any) 指定された引数の値、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
def getfloatarg(position, defval = None):
"""
概要:
コマンドライン引数を浮動小数点数に変換して取得します。
詳細説明:
getarg()で取得した文字列をpfloat()を使って浮動小数点数に変換します。
引数が存在しない、または浮動小数点数に変換できない場合、デフォルト値を返します。
:param position: (int) 取得する引数の位置（インデックス）。
:param defval: (float or None, optional) 引数が存在しないか、変換できない場合に返すデフォルト値。デフォルトはNone。
:returns: (float or None) 変換された浮動小数点数の引数、またはデフォルト値。
"""
return pfloat(getarg(position, defval))
def getintarg(position, defval = None):
"""
概要:
コマンドライン引数を整数に変換して取得します。
詳細説明:
getarg()で取得した文字列をpint()を使って整数に変換します。
引数が存在しない、または整数に変換できない場合、デフォルト値を返します。
:param position: (int) 取得する引数の位置（インデックス）。
:param defval: (int or None, optional) 引数が存在しないか、変換できない場合に返すデフォルト値。デフォルトはNone。
:returns: (int or None) 変換された整数の引数、またはデフォルト値。
"""
return pint(getarg(position, defval))
def usage(ka = ka, Z = Z, n = n, l = l, m = m):
"""
概要:
プログラムのコマンドライン引数の使用方法を表示します。
詳細説明:
現在のグローバル変数に基づいて、スクリプトの実行方法と利用可能なモードオプションを説明します。
:param ka: (float, optional) 軌道の指数係数。デフォルトはグローバル変数`ka`。
:param Z: (float, optional) 原子番号。デフォルトはグローバル変数`Z`。
:param n: (int, optional) 主量子数。デフォルトはグローバル変数`n`。
:param l: (int, optional) 方位量子数。デフォルトはグローバル変数`l`。
:param m: (int, optional) 磁気量子数。デフォルトはグローバル変数`m`。
:returns: 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()関数を呼び出して使用方法ガイドを表示し、
sys.exit()でプログラムを終了させます。
:param message: (str, optional) 表示するエラーメッセージ。デフォルトはNone。
:returns: None (プログラムが終了するため)
"""
if message is not None:
print("")
print(message)
usage()
print("")
exit()
def updatevars():
"""
概要:
コマンドライン引数に基づいてグローバル変数を更新します。
詳細説明:
sys.argvリストからコマンドライン引数を解析し、`mode`, `Z`, `ka`, `Ne` などの
グローバル変数を更新します。`mode`に'e'が含まれる場合、`ELabel`は'Etot'に設定され、
それ以外の場合は'E 1s'に設定されます。
:returns: 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) を計算します。
詳細説明:
規格化された形式 R_1s(r) = (2 / sqrt(4 * pi)) * (ka * Z)^1.5 * exp(-ka * Z * r)
を使用します。ここで、kaは指数部の係数、Zは有効原子核電荷、rは動径距離です。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数（この関数では1s軌道に対応するため1）。
:param l: (int) 方位量子数（この関数では1s軌道に対応するため0）。
:param m: (int) 磁気量子数（この関数では1s軌道に対応するため0）。
:param r: (float) 動径距離。
:returns: (float) 動径分布関数 R_nl(r) の値。
"""
return R1s0 * pow(ka * Z, 1.5) * exp(-ka * Z * r)
def rho(ka, Z, n, l, m, r):
"""
概要:
電子密度関数 rho(r) を計算します。
詳細説明:
動径分布関数 Rr(r) の絶対値の二乗として定義されます。
rho(r) = |Rr(r)|^2
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数（この関数では1s軌道に対応するため1）。
:param l: (int) 方位量子数（この関数では1s軌道に対応するため0）。
:param m: (int) 磁気量子数（この関数では1s軌道に対応するため0）。
:param r: (float) 動径距離。
:returns: (float) 電子密度関数 rho(r) の値。
"""
psi = Rr(ka, Z, n, l, m, r)
return psi * psi
def calQ(R):
"""
概要:
動径R内にある電荷量 Q(R) を計算します。
詳細説明:
事前にbuild_Qr()関数で計算・補間された電荷分布関数qfuncを用いて、
動径R内（0からRまで）に存在する全電子電荷量を返します。
Rが定義された範囲外の場合、エラーメッセージを出力しプログラムを終了します。
:param R: (float) 動径距離。
:returns: (float) 動径R内にある電荷量。
"""
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`における電子密度`rho(r)`を数値積分し、
各点での累積電荷量`Q(r)`を計算します。
計算された`Q(r)`のデータポイントを用いて、`scipy.interpolate.interp1d`により
3次スプライン補間を行い、グローバル変数`qfunc`として設定します。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:returns: None (グローバル変数`yRr`, `yQr`, `qfunc`を更新する)
"""
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-like 1s軌道の運動エネルギーの解析解 `Z^2 * ka^2 / 2 * (e^2 / (4 * pi * e0 * a0))` を計算します。
:param Z: (float, optional) 有効原子核電荷。デフォルトは1.0。
:param ka: (float, optional) 軌道の指数係数。デフォルトは1.0。
:returns: (float) 計算された運動エネルギー（Hartree単位）。
"""
# 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-like 1s軌道の核-電子間ポテンシャルエネルギーの解析解 `-Z^2 * ka^2 * (e^2 / (4 * pi * e0 * a0))` を計算します。
:param Z: (float, optional) 有効原子核電荷。デフォルトは1.0。
:param ka: (float, optional) 軌道の指数係数。デフォルトは1.0。
:returns: (float) 計算された核-電子間ポテンシャルエネルギー（Hartree単位）。
"""
# 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):
"""
概要:
台形公式を用いて関数を数値積分します（現在未使用）。
詳細説明:
与えられた関数をE0からE1まで、ステップサイズhで台形公式を使用して数値積分します。
この関数は現在メインロジックでは使用されていません。
:param func: (callable) 積分対象の関数。
:param E0: (float) 積分範囲の下限。
:param E1: (float) 積分範囲の上限。
:param h: (float) ステップサイズ。
:returns: (list) [積分結果, -1.0] (エラー値は現在ダミーの-1.0)。
"""
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次元の動径積分 integral(4pi * r^2 * func(r) dr) を実行します。
詳細説明:
球座標における動径成分の積分を`scipy.integrate.quad`を用いて計算します。
積分対象の関数`func`は`f(r)`の形式で与えられます。
:param func: (callable) 積分対象の関数 f(r)。
:param rmin: (float) 積分範囲の下限。
:param rmax: (float) 積分範囲の上限。
:param limit: (int, optional) `quad`関数のサブインターバル数の最大値。デフォルトは15。
:param epsrel: (float, optional) 積分結果の相対誤差。デフォルトは1.0e-8。
:returns: (tuple) (積分結果 I, 推定誤差 err)。
"""
I, err = integrate.quad(lambda r: pi4 * r * r * func(r), rmin, rmax, limit, epsrel)
return I, err
def calUZ(r, Z):
"""
概要:
原子核による静電ポテンシャル U_Z(r) を計算します。
詳細説明:
点電荷である原子核（電荷 Z）が動径rの位置に作るクーロンポテンシャル `-Z / r` を計算します。
rが非常に小さい場合（1.0e-3未満）は、数値的な安定性のためrを1.0e-3として扱います。
:param r: (float) 動径距離。
:param Z: (float) 原子番号（核電荷）。
:returns: (float) 核による静電ポテンシャル。
"""
if r < 1.0e-3:
return -Z / 1.0e-3
else:
return -Z / r
def calUrho(r, ka, Z, n, l, m):
"""
概要:
電子密度 rho(r) によって形成される静電ポテンシャルを計算します。
詳細説明:
動径rの位置における電子-電子間の静電ポテンシャルを計算します。
このポテンシャルは、動径r内にある電荷Q(r)によるポテンシャルQ(r)/rと、
動径r外にある電荷によるポテンシャル積分項の和として構成されます。
:param r: (float) ポテンシャルを計算する動径距離。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:returns: (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):
"""
概要:
運動エネルギーを計算します。
詳細説明:
1s軌道の動径関数 Rr(r) を用いて、ラプラシアン演算子を含む運動エネルギーの積分を数値的に計算します。
積分は解析的に簡略化された形 `(1/2) * 4pi * Z * ka * R1s0^2 * (ka * Z)^3 * integral(r * (2 - Z * ka * r) * exp(-2 * Z * ka * r) dr)`
として実行されます。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数（この関数では直接使用されない）。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (運動エネルギー T, 積分誤差 errT)。
"""
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):
"""
概要:
核-電子間ポテンシャルエネルギーを計算します。
詳細説明:
電子密度 rho(r) と原子核のクーロンポテンシャル (-Z/r) との相互作用エネルギーを計算します。
積分は `integral(4pi * r^2 * (-Z/r) * rho(r) dr)` の形式で実行されます。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (核-電子間ポテンシャルエネルギー UeZ, 積分誤差 errUeZ)。
"""
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):
"""
概要:
電子-電子間反発エネルギーを計算します。
詳細説明:
電子密度 rho(r) によって形成される静電ポテンシャル Urho(r) と、
電子密度 rho(r) との相互作用エネルギーを計算します。
積分は `integral(4pi * r^2 * Urho(r) * rho(r) dr)` の形式で実行されます。
`mode`に'e'が含まれる場合、全エネルギー計算モードとみなされ、Neの二乗が係数として適用されます。
それ以外の場合は、1電子エネルギーモードとみなされ、Neが係数として適用されます。
:param mode: (str) 計算モード。'e'が含まれる場合、全エネルギーモードとしてNe*Neが乗算される。
それ以外の場合はNeが乗算される。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (電子-電子間反発エネルギー Uee, 積分誤差 errUee)。
"""
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-alpha交換ポテンシャルの局所値を計算します。
詳細説明:
動径rにおけるX-alpha交換ポテンシャル `(-3.0 * alpha * (3 / (4 * pi) * rho(r))^(1/3))` を計算します。
ここで`alpha`はSlaterの交換相関係数です。
:param r: (float) 動径距離。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数（この関数では直接使用されない）。
:returns: (float) X-alpha交換ポテンシャルの局所値。
"""
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-alpha交換エネルギーを計算します。
詳細説明:
X-alpha交換ポテンシャル `calLocalUXa(r)` と電子密度 `rho(r)` との相互作用エネルギーを計算します。
積分は `integral(4pi * r^2 * calLocalUXa(r) * rho(r) dr)` の形式で実行されます。
`mode`に'e'が含まれる場合、全エネルギー計算モードとみなされ、Neの4/3乗が係数として適用されます。
それ以外の場合は、1電子エネルギーモードとみなされ、Neの1/3乗が係数として適用されます。
:param mode: (str) 計算モード。'e'が含まれる場合、全エネルギーモードとしてNe^(4/3)が乗算される。
それ以外の場合はNe^(1/3)が乗算される。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (X-alpha交換エネルギー UXa, 積分誤差 errUXa)。
"""
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):
"""
概要:
全エネルギー Etot を計算します。
詳細説明:
運動エネルギー (T)、核-電子間ポテンシャルエネルギー (UeZ)、電子-電子間反発エネルギー (Uee)、
およびX-alpha交換エネルギー (UXa) の各成分を計算し、その合計として全エネルギーを算出します。
計算前に`build_Qr`を呼び出して電荷分布関数を構築します。
:param ka: (float) 1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (T, UeZ, Uee, UXa, Etot) - それぞれのエネルギー値と全エネルギー。
"""
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):
"""
概要:
全エネルギー Etot の値のみを計算して返します。
詳細説明:
calTotalEnergyと同様に各エネルギー成分を計算し合計しますが、
主に`scipy.optimize.minimize`などの最適化関数に渡すために、
単一の全エネルギー値を返します。この関数内では`mode`は'e'として扱われます。
:param kav: (float) 1s動径関数の指数部の係数（最適化される変数）。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (float) 計算された全エネルギー Etot。
"""
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`を用いて、`calTotalEnergyOnly`関数を目的関数として、
`ka`を最適化します。選択された最適化手法（例: Nelder-Mead、CGなど）に基づいて、
最小エネルギーを与える`ka`の値を見つけ出します。
最適化された`ka`を用いて、最終的なエネルギー成分と全エネルギーを再計算して返します。
:param ka: (float) 初期値としての1s動径関数の指数部の係数。
:param Z: (float) 有効原子核電荷。
:param n: (int) 主量子数。
:param l: (int) 方位量子数。
:param m: (int) 磁気量子数。
:param Ne: (float) 電子数。
:param Rmin: (float) 積分範囲の最小値。
:param Rmax: (float) 積分範囲の最大値。
:param nmaxdiv: (int) 数値積分の分割数。
:param eps: (float) 数値積分の相対誤差。
:returns: (tuple) (T, UeZ, Uee, UXa, Etot, optimized_ka) - 各エネルギー成分、全エネルギー、および最適化された`ka`。
"""
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):
"""
概要:
最適化アルゴリズムで使用する一次導関数を前方差分法で近似計算します。
詳細説明:
`scipy.optimize.minimize`で勾配情報が必要なアルゴリズムのために、
`calTotalEnergyOnly`を呼び出して、パラメータ`ai`の各成分に対して
前方差分近似を用いて数値的に勾配を計算します。
差分幅はグローバル変数`h_diff`で設定されます。
:param ai: (numpy.ndarray) 最適化されるパラメータの配列（ここでは`ka`のみを含む）。
:returns: (numpy.ndarray) 各パラメータに対する一次導関数（勾配）の配列。
"""
global Z, n, l, m, Ne, Rmin, Rmax, nmaxdiv, eps
n = len(ai)
f0 = Emin = 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() のコールバック関数です。
詳細説明:
最適化の各イテレーションで呼び出され、現在のパラメータ値 (xk) と
それに対応する全エネルギーをコンソールに出力します。
これは最適化の進行状況を監視するために使用されます。
:param xk: (numpy.ndarray) 現在の最適化パラメータの配列。
:returns: 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 を掃引し、各Neにおけるエネルギー成分と全エネルギーを計算・出力します。
詳細説明:
グローバル変数`Nearray`に定義された電子数`Ne`の各値について、
全エネルギーおよび運動エネルギー、核-電子間、電子-電子間、X-alpha交換の各エネルギー成分を計算します。
計算結果は表形式で出力され、オプションで`ka`を最適化した結果も計算・出力し、
さらにエネルギー対Neのグラフをプロットします。
Ne=0.5付近でのエネルギーの放物線展開も計算・表示されます。
:param Eanal: (float) 解析的に導出された全エネルギー（比較用）。
:returns: 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 を掃引し、各kaにおけるエネルギー成分と全エネルギーを計算・出力します。
詳細説明:
グローバル変数`kaarray`に定義された軌道指数係数`ka`の各値について、
全エネルギーおよび運動エネルギー、核-電子間、電子-電子間、X-alpha交換の各エネルギー成分を計算します。
計算結果は表形式で出力され、オプションでエネルギー対kaのグラフをプロットします。
:param Eanal: (float) 解析的に導出された全エネルギー（比較用）。
:returns: None
"""
global Rmin, Rstep, nR, Rmaxdata, r
global qfunc
global nmaxdic, 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):
"""
概要:
デバッグモードで単一の計算結果と各種関数・ポテンシャルをプロットします。
詳細説明:
現在のグローバル変数（ka, Z, Neなど）に基づいて一度全エネルギー計算を実行し、
その結果（各エネルギー成分と全エネルギー）を表形式で出力します。
モードに'v'が含まれる場合、kaを最適化した結果も同様に出力します。
さらに、動径分布関数 Rr(r)、動径R内電荷 Q(r)、核ポテンシャル U(Z)、
核-電子ポテンシャル U(Z-rho)、電子-電子ポテンシャル -U(rho)、X-alphaポテンシャル U(Xa)
のグラフをプロットします。
:param Eanal: (float) 解析的に導出された全エネルギー（比較用）。
:returns: 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():
"""
概要:
プログラムのメインエントリポイントです。
詳細説明:
コマンドライン引数からグローバル変数を更新し、動径グリッドを初期化します。
H-like原子の解析解を計算・表示し、数値積分の動径関数規格化チェックを行います。
その後、`mode`変数に応じて、デバッグモード(`d`)、ka掃引モード(`k`)、
またはNe掃引モード(`n`)のいずれかの処理を実行します。
:returns: 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("epsR=", epsR, Z, ka)
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()