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

簡単のため、二個の離れた原子核からのクーロンポテンシャルのみを扱うとする。
この時の、厳密なポテンシャルと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()