簡単のため、二個の離れた原子核からのクーロンポテンシャルのみを扱うとする。
この時の、厳密なポテンシャルとMuffin-tin近似との差を見る。
Muffin-tin近似をする際、隣のサイトのポテンシャルの球平均は、以前にまとめた方法を使った。
koideforest.hatenadiary.com
import numpy as np from matplotlib import pyplot as plt def norm( x, y, z, x0, y0, z0 ): vec = [ x - x0, y - y0, z - z0 ] return np.linalg.norm( vec ) def step( x, y, z, x0, y0, z0, rmt ): if norm( x, y, z, x0, y0, z0 ) > rmt: return 0. else: return 1. def site_v( x, y, z, x0, y0, z0 ): return - 1. / ( norm( x, y, z, x0, y0, z0 ) + 1e-7 ) def V( x, y, z ): xa = 2. return site_v( x, y, z, 0, 0, 0 ) + site_v( x, y, z, xa, 0, 0 ) def site_vmt( x, y, z, x0, y0, z0 ): R = 2. pot_next = - 1. / R pot = - 1. / ( norm( x, y, z, x0, y0, z0 ) + 1e-7 ) + pot_next return pot def VMT( x, y, z ): xa = 2. R = xa step1 = step( x, y, z, 0, 0, 0, R/2 ) step2 = step( x, y, z, xa, 0, 0, R/2 ) VMT0 = site_vmt( xa/2, 0, 0, 0, 0, 0 ) * ( 1 - step1 ) * ( 1 - step2 ) pot1 = site_vmt( x, y, z, 0, 0, 0 ) * step1 pot2 = site_vmt( x, y, z, xa, 0, 0 ) * step2 return pot1 + pot2 + VMT0 N = 1000 x = np.linspace( -5, 5, N ) y = np.linspace( -5, 5, N ) X, Y = np.meshgrid( x, y ) mat_v = np.zeros( ( N, N ) ) mat_vmt = np.zeros( ( N, N ) ) for ix, xx in enumerate( x ): for iy, yy in enumerate( y ): mat_v[ ix ][ iy ] = V( xx, yy, 0 ) mat_vmt[ ix ][ iy ] = VMT( xx, yy, 0 ) fig = plt.figure( figsize = ( 5, 4 ) ) ax = fig.add_subplot( 1, 1, 1 ) pcm = ax.pcolormesh( X, Y, mat_v.T, vmin = -3, vmax = 0 ) fig.colorbar( pcm ) plt.show() fig = plt.figure( figsize = ( 5, 4 ) ) ax = fig.add_subplot( 1, 1, 1 ) pcm = ax.pcolormesh( X, Y, mat_vmt.T, vmin = -3, vmax = 0 ) fig.colorbar( pcm ) plt.show() fig = plt.figure( figsize = ( 5, 4 ) ) ax = fig.add_subplot( 1, 1, 1 ) pcm = ax.pcolormesh( X, Y, mat_v.T - mat_vmt.T, cmap='bwr', vmin = -1.5, vmax = 1.5 ) fig.colorbar( pcm ) plt.show()
結合方向のある狭い範囲で、ポテンシャルの差が大きいことがわかる。
等高線プロットで参考にしたページ。
[Pythonによる科学・技術計算] 2次元(カラー)等高線等の描画,可視化,matplotlib - Qiita