浜本さんのコード

コードを理解／移植する場合のメモ
- おそらく計算時間を短縮に重点が置かれているようだ。そのために、いろいろテクニックを使っている。
- ポテンシャルの近くでは、計算点を多くとり計算精度をあげている。r = 0.025, 0.050, 0.075, 0.100, 0.125, 0.150, 0.175, 0.200, 0.250, 0.300, 0.350, 0.400, 0.500, 0.600, 0.700, 0.800, 1.000, 1.2000, 1.400, ... fm.
- 実数変数は Double precision を用いているが、実数値は単精度。例えば 1.0/EX など。すべて倍精度にするには、1.0d0/EX と書き換える。
- ファイルへの出力は小数点以下6桁で出力しているので、場合によっては（コードの移植のために計算値をプロットして比較するときなど）精度が足りない。
- 59.4 fm からおよそ 3.8 fm まで 0.2 fm 間隔で Störmer method を用いてシュレディンガー方程式を積分していく(参考: Encyclopedia of Mathematics - Störmer method)。積分の精度をあげると (例えば 0.2 fm を 0.1 fm にする) V_WS の値が 100 keV 程度変わる。³¹Ne の場合、S_n = -300keV の p_3/2 軌道を計算すると、0.2 fm の精度の場合 V_WS = -46.75 MeV となるが (PRC81(2010)021304)、精度をあげていくと V_WS = -46.83 MeV あたりに近づくようだ。
- Störmer method は前 5 点の二階微分の値が必要。コード内だと r = 0.0 fm における値を 4 回使っている。より精度の高い計算をすると、この初期値が効いてくる。
- 外側から計算した波動関数と内側から計算した波動関数を 3.8 fm 地点で接続させる際、傾きを計算するが、前後 2 点, 計 4 点の値を使って正確に計算しているようだ。しかし、ポテンシャルの深さの決定には、ここはあまり効かないかも？
- コードを移植する際、同じアルゴリズムで倍精度の変数を用いれば、大体の値は 9 桁以上で一致する。ただし、シュレディンガー方程式を端から積分していく場合、最初は9桁で合っていた値でも若干変わっていき、5桁程度の一致になる（r = 59.8 fm --> 3.8 fm の場合）。
main-eg.f

サブルーチン	変数名	意味	値
MAIN	ISTRTW	バレンス中性子のスピンの大きさを 2 倍したもの	1
MAIN	VSO	dvipng error! exitcode was 2 (signal 0), transscript follows: の v　(PRC69(2004)041306(R). (1) 式参照) 。 V_so(r) の方ではないので注意。	32 (Β安定線上の核のスタンダードな値)
MAIN	DFN	diffuseness。 dvipng error! exitcode was 2 (signal 0), transscript follows: の a (PRC69(2004)041306(R). (2) 式参照)。	0.67 fm (Β安定線上の核のスタンダードな値)
MAIN	RZERO	dvipng error! exitcode was 2 (signal 0), transscript follows: の r₀ (PRC69(2004)041306(R). (2) 式参照)。	1.27 fm (Β安定線上の核のスタンダードな値)
MAIN	TMAS	コアの質量数。これを A_T とする。	-
MAIN	PMAS	バレンス中性子の質量数。これを A_P とする。	1
MAIN	RMAS	コアとバレンス中性子の質量数の換算質量"数"。これを A_μ とする。 1/A_μ = 1/A_P + 1/A_T → A_μ = A_TA_P/(A_T+A_P)	大雑把には 1
MAIN	CHARGE	バレンス中性子の電荷 Z	0
MAIN	WNUNIT	dvipng error! exitcode was 2 (signal 0), transscript follows: 。M は中性子または陽子の質量 (変数 CHARGE によって変わる)。	中性子 : 0.21968 MeV^-1/2fm^-1 陽子 : 0.21953 MeV^-1/2^fm-1^
MAIN	CFUNIT	dvipng error! exitcode was 2 (signal 0), transscript follows: 。M は中性子または陽子の質量。	大雑把には ?
MAIN	XBFAC	A_T^1/3	-
MAIN	XBAR	R = r₀A_T^1/3。コアの核半径。	-
MAIN	XCPLE	動径方向の波動関数を出力するときの範囲	60 fm
MAIN	XMES2	動径方向の波動関数を出力するときのステップ。これを Δr とする。	0.2 fm
MAIN	NXMAX	計算範囲のステップ数	300
MAIN	NXRM	コアの核半径のステップ数。 R/Δr。	-
MAIN	RCWL	dvipng error! exitcode was 2 (signal 0), transscript follows: 。Reduced Compton WaveLength of nucleon dvipng error! exitcode was 2 (signal 0), transscript follows: を 2 で割ったもの dvipng error! exitcode was 2 (signal 0), transscript follows: 。	-
MAIN	JJ1, JJROW, J1	バレンス中性子の全角運動量 j の大きさを 2 倍したもの	-
MAIN	LL1, LLROW, L1	バレンス中性子の軌道角運動量 l の大きさを 2 倍したもの	-
VMACH	SPFC	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
VMACH	SPFACT	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
VMACH	EZERO, ECM	中性子分離エネルギー S_n のマイナス	-
VMACH	VSX	ウッズサクソンポテンシャルの深さ V_WS の初期値	典型的には -40~50 MeV
VMACH	DELEDR	ウッズサクソンポテンシャルの深さ V_WS をサーチする際の微小変化量	0.0001 MeV
POTENV	COE1	サブルーチン POTENV 内で用いるウッズサクソンポテンシャルの深さ V_WS	-
POTENV	COE2	Z_eff = 1.4398650 Z (?)	-
POTENV	COE3	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	COE4	V_WS/a	-
POTENV	XBARC	R_c = r_0''c''A_T^1/3。コアの陽子半径。	-
POTENV	VCLFC1	Z_eff/(2*R_c)	-
POTENV	NDFMES	原点付近の波動関数を精度良く計算するため、動径方向を分割するが、その時の分割領域数。	4
POTENV	XMES1	原点付近での波動関数を計算する際のステップ	-
POTENV	X	動径方向の変数 r	-
POTENV	DX	動径方向の変数の微小区間 Δr	-
POTENV	EX	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	EXINV	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	NX	各微小区間の番号	1 ~ 300?
POTENV	VCENTR(NX)	各微小区間に対する dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	VSPIN(NX)	各微小区間に対する dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	VCOULM(NX)	各微小区間に対するクーロンポテンシャル。 dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	VKR(NX)	各微小区間に対する dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	VKRS1(NX)	各微小区間に対する dvipng error! exitcode was 2 (signal 0), transscript follows:	-
POTENV	VKRS2(NX)	各微小区間に対する dvipng error! exitcode was 2 (signal 0), transscript follows:	-
VMACH	E1	中性子分離エネルギー S_n の絶対値	-
VMACH	WN	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
VMACH	WNINI	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
AHNKEV	L1	バレンス中性子の軌道角運動量 l の大きさを 2 倍したもの	-
AHNKEV	LL	バレンス中性子の軌道角運動量 l の大きさ + 1	-
AHNKEV	X	計算範囲の少し外 (NXMAX+3) * XMES2 = 60.6 からスタート。	-
AHNKEV	NQ	((NXMAX+2)XMES2 ~ (NXMAX-2)XMES2) : 60.4, 60.2, 60.0, 59.8, 59.6 fm に対応する 5 点に対して Hankel 関数を計算。	-
AHNKEV	RHOMX,R	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
AHNKEV	RI	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
AHNKEV	EX	dvipng error! exitcode was 2 (signal 0), transscript follows:	-
AHNKEV	G(1)	EX * RI = dvipng error! exitcode was 2 (signal 0), transscript follows:	-
AHNKEV	G(2)	EX * (RI + RI**2)	-
AHNKEV	G(3)	EX * (2RI + 3RI*2 + 3RI**3)	-
AHNKEV	G(L+2)	G(L+2) = (2L + 1) * RI * G(L+1) + G(L)	-
AHNKEV	LL2	バレンス中性子の軌道角運動量 l の大きさ - 1	-
AHNKEV	GCRA(LL,NQ)	バレンス中性子の軌道角運動量 l の大きさ + 1 と 60.4, 60.2, 60.0, 59.8, 59.6 fm の 5 点に対する Hankel 関数の値を格納？	-

サブルーチン VMACH での各変数の値

ITEI	INEX	NDMAX	ND	FND	DX	K2COR2	X_ini	NXIMIN	N3	X	NXIMAX	N3	X	NXMAX3	NXPRCH	NNX
1	1	4	1	1	0.025	-1	0	2	1	0.025	8	7	0.175		3	2
1	1	4	2	2	0.05	3	0.15	5	8	0.2	8	11	0.35		5	3
1	1	4	3	4	0.1	7	0.3	5	12	0.4	8	15	0.7		5	3
1	1	4	4	8	0.2	11	0.6	5	16	0.8	21	32	4	17	5	3
1	2	1	1	1	-0.2	311	59.8	1	310	59.6	281	30	3.6	276

動径方向シュレディンガー方程式の解き方

dvipng error! exitcode was 2 (signal 0), transscript follows:

Störmer method を用いて解を求めている(参考: Encyclopedia of Mathematics - Störmer method)。これで解くのはなにか理由があるのだろうか？

cc-nilson-g.f

dvipng error! exitcode was 2 (signal 0), transscript follows:

ここで

dvipng error! exitcode was 2 (signal 0), transscript follows:

より(これは違ってた。。。以下間違ってるかも)、
- ```
dvipng error! exitcode was 2 (signal 0), transscript follows:
```

hamamoto_code

おぼえがき

リンク

浜本さんのコード