nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

プログラミング

FFTWで3次元高速フーリエ変換してみた。

前回の例を元に、3次元に拡張してみた。 koideforest.hatenadiary.com program fftw3_test_3d implicit none integer, parameter :: WP = selected_real_kind( 8 ) real( WP ), parameter :: PI = 3.141592653589793 complex( WP ), parameter :: I = ( 0, 1…

FFTWをfortranで使う際の自分用のメモ

FFTWをfortranで使う。本家の説明: FFTW 3.3.8: Fortran Examples参考にしたサイト: fftwの使い方:腰も砕けよ 膝も折れよ:So-net blog FortranのFFTWの使い方 - Qiita Hiroyuki R. Takahashi Fumio KUSUHARA -Fortran90で構造解析- program fftw3_test i…

Fortranの精度変更。

Fortranの精度変更の書式が自分の知らない感じだったので、自分用にまとめる。 Fortran Best Practices — Fortran90 1.0 documentation program test_presicion implicit none integer, parameter :: dp0 = kind( 0.d0 ) integer, parameter :: dp = selecte…

Makefileの自分用のメモ

FC = gfortran FFLAGS = INCLUDE = -I/opt/local/include LIBS = -L/opt/local/lib -lbrabra PROGRAM = main.out OBJS = sub.o SRCS = $(OBJS:%.o=%.f90) .Phony: all clean all: $(PROGRAM) $(PROGRAM): $(OBJS) $(FC) $(FFLAGS) -o $@ $^ $(INCLUDE) $(LIB…

pythonで三次元プロットする時の自分用メモ

三次元プロット用のx軸、y軸の配列。 import numpy as np N = 100 x_min, x_max = -1, 1 y_min, y_max = -1, 1 x1 = np.linspace( x_min, x_max, N ) y1 = np.linspace( y_min, y_max, N ) X, Y = np.meshgrid( x1, y1 ) 三次元プロット用のz軸データ。 def …

自分が自分じゃないものってNaNだ。

"NaN"とは"Not a Number"の略で、例えば「1/0」のような変な計算をしてしまった時に出力される。 このNaNの特徴として、「NaN = NaN」が「偽(False)」になる。 逆に言うと、「NaN is not NaN」が「真(True)」になる。NaNとは「自分が自分でないもの」な…

波動関数の節を数える

文字通り、波動関数の節を数える。ndarrayの中で条件を満たす要素数を数える方法。 NumPy配列ndarrayの条件を満たす要素数をカウント | note.nkmk.me import numpy as np x_min, x_max, N = 0, 3*np.pi, 100 x = np.linspace( x_min, x_max, N ) wave = np.c…

Spigot Algorithm

ネイピア数ってそういえばどうやって計算しているのかと思って調べてみたら、面白い記事を見つけた。 こつこつアルゴリズム(Spigot Algorithm)無限級数が以下のように定義されているとする。 この時、 と書けるので、 ここで、 したがって、これを次々繰り返…

二階微分演算子の変数変換 for Numerov method

Numerov法は一階微分方程式を含まない二階微分方程式を解くのに便利な方法である。 この方法を適用するには、等間隔の変数刻み(メッシュ or グリッド)が必要である。 しかし、変化が穏やかなところでは粗いメッシュで十分だが、一方で変化が急なところでは…

matplotlibで二文字以上の変数を上付きにする。

matplotlibで普通に from matplotlib import pyplot as plt atomic_number = 29 plt.title( 'Cu$^{}$'.format( atomic_number ) ) とやると、「Cu9」となってしまう。これは、内部で「$^29$」と書かれたと見做されているためである。 ベタ打ちするならば、「…

ポテンシャルのルジャンドル関数展開

「ポテンシャルの角運動量展開」の二次元版だと思って頂いて差し支えない。係数 はルジャンドル関数の直交関係から出て来る。 koideforest.hatenadiary.com上記の式を用いて、まで展開してみる。 import numpy as np from math import radians from scipy.sp…

Muffin-tin近似と厳密なポテンシャルとの差について

簡単のため、二個の離れた原子核からのクーロンポテンシャルのみを扱うとする。 この時の、厳密なポテンシャルとMuffin-tin近似との差を見る。Muffin-tin近似をする際、隣のサイトのポテンシャルの球平均は、以前にまとめた方法を使った。 koideforest.haten…

誤差関数によるステップ関数でGibbs現象は起こるか?

ステップ関数等で不連続に打ち切られた関数をフーリエ変換しようとすると、どんなに頑張っても振動が残る。 これはギブス現象として知られている。 ギブズ現象 - Wikipediaでは、ステップ関数の代わりに誤差関数で滑らかにしたら、どれくらい収束が良いのか…

Juliaで一次元井戸型ポテンシャル

以下のサイトの下の方に、Juliaで一次元のシュレーディンガー方程式を解くPDFが紹介されている。 物理ノートby永井Juliaの練習としてやってみた。 PDF内では、無限の井戸の中に斥力ポテンシャルを入れた場合をやっているが、ここでは引力ポテンシャルに対し…

pythonでペンローズタイリング

ペンローズタイリングとは、ペンローズが開発した非周期的な図形の敷き詰め方である。 以下のサイトで、ペンローズタイリングの方法と、そのpythonスクリプトが公開されている。 Penrose Tiling Explained図の描画には、pycairoが使われていて、事前にインス…

頂角が小さい二等辺三角形の底辺について

回転運動の説明を読んでいる時に、角度が小さい時の差ベクトルの近似について気になったので考えて見た。 等速円運動 [物理のかぎしっぽ]以下の図のような二等辺三角形における底辺(青)、円弧(赤)、そして垂線(緑)を考える。 図では頂角は30度でプロッ…

sympyで(平面の)Greenの定理を確認

sympyの練習を兼ねて、平面に対するGreenの定理で、問題を解いてみる。 平面のグリーンの定理 [物理のかぎしっぽ]以下の積分を求めてみる。 原点で最大値を取り、等方向的で、境界でゼロを持つような、何かしらの密度を積分する、と思えば解りやすいだろうか…

多重項を数える(2)

以前に、軌道内の多重項を求めるスクリプトを書いた。 koideforest.hatenadiary.comここではそれらを用いた上で、さらに軌道間で多重項を合成する。 軌道の記号と値を行き来する関数 def spec2l( spec ): if spec in { "s", "S" }: l = 0 elif spec in { "p"…

多重項を数える(1)

前回、binary表現を使って、電子配置を作ることを試みた。 koideforest.hatenadiary.com次のステップとして多重項を数え上げるため、電子配置を作るのと同時に角運動量とスピンも一緒にリストで保存するように改良した。 def add_electron( configurations0,…

バイナリーと量子力学

例えば2p軌道とかのp軌道の場合、一電子状態は以下のどれかに対応する。(1,u), (1, d), (0, u), (0, d), ( -1, u), ( -1, d)これを、binary形式、「電子がいる軌道=1」、「電子がいない軌道=0」とすれば、例えば、(1,u)と(1, d)が占有されているとき、110000…

XMCDのSum ruleをpythonで処理

とあるスクールで、講師が事前に用意したOriginのマクロを使ってXMCDのSum ruleから磁気モーメントを求めるという講習があったが、Originを使いたいと全く思わなかったため、一人でpythonスクリプトを黙々と作っていた。ローレンツ関数やら誤差関数やらで適…

MacOSでpythonをインストールし直す

python2を主に使っていたが、仕事上ちゃんとpython3の環境を整えなければならなくなり、その際の記録。いろいろお試し感覚でごちゃごちゃ入れていて、 System由来のpython2 手動経由のpython2 MacPorts経由のpython2 pyenv経由のmicroconda (python3) が混在…

横軸と縦軸のデータをすぐに保存&読込 (python)

pythonで作業していて、パッとデータを保存する、もしくは読み込む方法。 結局は、横軸と縦軸の各一次元データ(横ベクトル)を結合して縦長のデータにする部分が面倒くさい。 しかし、それは numpy.c_ で解決することがわかった。 Pythonで配列や行列の結合…

ガウス関数のフーリエ変換 (scipy.fftpack)

ガウス関数はフーリエ変換しても、数式上は(幅は変わるが)またガウス関数になる。 高速フーリエ変換(FFT)とフーリエ変換との関係を以前まとめたため、FFTを用いてガウス関数をフーリエ変換する。 koideforest.hatenadiary.comここではscipy.fftpackを用…

散乱振幅からサイトT行列を復元する。

実験結果か何かで、散乱振幅の角度分布が得られているとする。 ここからサイト行列を復元する方法を考える。 サイト行列が求まれば、そこからArgand diagramを作って、(準)共鳴準位があるとかないとか、追加で情報が得られる(もし角度分布のエネルギー依…

部分散乱振幅(サイトT行列)のアーガンド図

仕事で図を作る必要があったので、練習がてらまとめた。位相シフトなどの概念にについては、wikipedia等を参照して頂きたい。 散乱振幅 - Wikipedia 簡単に言うと、「ポテンシャルが無い時の波は、ポテンシャルによってどれくらい変調されたか?」を表す量で…

水素様原子波動関数 (python)

以前、pythonにおける水素様原子のためのLaguerre陪関数について確認した。 koideforest.hatenadiary.comせっかくなので、pythonで規格化まで含めた水素様原子波動関数のスクリプト(原子単位系)を作った。 from scipy.special import assoc_laguerre from …

X線吸収の多重散乱理論でよく使いそうな関数 (python)

ほぼ自分用のメモ from scipy.special import spherical_jn, spherical_yn, sph_harm from sympy import * from sympy.physics.wigner import clebsch_gordan, gaunt # spherical Bessel function def jl( l, rho ): return spherical_jn( l, rho ) # spheri…

EXAFS振動のフーリエ変換について

「EXAFS振動をフーリエ変換する」というフレーズはX線吸収分光をやっていると腐る程見かけるが、実際にはArtemisやらLarchやらのソフトで流し込んで、実際に気にするところはパラメータのフィッティングだけということがほとんどな気がする。 ここでは「そも…

メタン(CH4)の分子軌道とLiebフラットバンド

四つ足で御馴染みのメタンの分子軌道が、固体におけるLiebフラットバンドが現れる構造と同じということで、テンションが上がった。以下の簡単な飛び移りだけを考えたハミルトニアン行列の固有値と固有ベクトルを求める。 対角成分(各原子上でのエネルギー)…