crystal_MP_simple.py ダウンロード/コピー

crystal_MP_simple.py をダウンロード

crystal_MP_simple.py
crystal_MP_simple.py
  1"""
  2概要: マデルングポテンシャルを単純総和法で計算するスクリプト。
  3詳細説明:
  4    このスクリプトは、指定された結晶構造（格子定数とサイト情報）に基づき、
  5    単純総和法を用いてマデルングポテンシャルを計算します。
  6    計算されたポテンシャルは、中心イオンからの距離の関数としてプロットされます。
  7    tkcrystalbase.pyモジュールに依存します。
  8
  9関連リンク:
 10    crystal_MP_simple_usage
 11"""
 12import sys
 13import os
 14from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
 15import numpy as np
 16from numpy import linalg as la
 17import matplotlib.pyplot as plt
 18
 19from tkcrystalbase import *
 20
 21
 22pi          = 3.14159265358979323846
 23pi2         = pi + pi
 24torad       = 0.01745329251944 # rad/deg";
 25todeg       = 57.29577951472   # deg/rad";
 26basee       = 2.71828183
 27
 28h           = 6.6260755e-34    # Js";
 29h_bar       = 1.05459e-34      # "Js";
 30hbar        = h_bar
 31c           = 2.99792458e8     # m/s";
 32e           = 1.60218e-19      # C";
 33me          = 9.1093897e-31    # kg";
 34mp          = 1.6726231e-27    # kg";
 35mn          = 1.67495e-27      # kg";
 36u0          = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
 37e0          = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
 38e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
 39a0          = 5.29177e-11      # m";
 40kB          = 1.380658e-23     # JK<sup>-1</sup>";
 41NA          = 6.0221367e23     # mol<sup>-1</sup>";
 42R           = 8.31451          # J/K/mol";
 43F           = 96485.3          # C/mol";
 44g           = 9.81             # m/s2";
 45
 46
 47
 48# Lattice parameters (angstrom and degree)
 49#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
 50lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
 51
 52# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
 53sites = [
 54         ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 55        ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.5, 0.5])]
 56        ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.0, 0.5])]
 57        ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.0])]
 58        ,['Cl', 'Cl1', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
 59        ,['Cl', 'Cl2', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
 60        ,['Cl', 'Cl3', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
 61        ,['Cl', 'Cl4', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
 62        ]
 63
 64# r range
 65rmin =   0.1
 66rmax = 100.0
 67nr   = 101
 68
 69
 70# Figure configuration
 71figsize = (6, 6)
 72
 73rstep = (rmax - rmin) / (nr - 1)
 74
 75
 76def usage():
 77    """
 78    概要: スクリプトの正しい使用方法を表示します。
 79    詳細説明: コマンドライン引数rmaxとnrの指定方法をユーザーに示します。
 80    :returns: なし
 81    :rtype: None
 82    """
 83    print("")
 84    print("Usage: python {} rmax nr".format(argv[0]))
 85    print("   ex: python {} {} {}".format(argv[0], rmax, nr))
 86    print("")
 87
 88def terminate():
 89    """
 90    概要: エラー発生時に使用方法を表示してスクリプトを終了します。
 91    詳細説明: usage関数を呼び出して使用方法を出力した後、システムを終了します。
 92    :returns: なし (スクリプトが終了するため)
 93    :rtype: None
 94    """
 95    usage()
 96    exit()
 97
 98
 99def draw_box(ax, aij, nrange, color = 'black'):
100    """
101    概要: 結晶の単位格子境界を3Dプロットに描画します。
102    詳細説明: 単位格子の辺を黒線（または指定された色）で描画します。
103              この関数はdraw_unitcellから呼び出されますが、nrangeは現在の実装では使用されていません。
104    引数:
105        :param ax: matplotlibの3D軸オブジェクト。
106        :type ax: matplotlib.axes._subplots.Axes3DSubplot
107        :param aij: (3, 3)のndarray、格子ベクトルa, b, cを表す。
108        :type aij: numpy.ndarray
109        :param nrange: 描画する単位格子の範囲。[[xmin, xmax], [ymin, ymax], [zmin, zmax]]の形式。（現状未使用）
110        :type nrange: list of list of float
111        :param color: 描画する線の色。デフォルトは'black'。
112        :type color: str
113    戻り値:
114        :returns: なし
115        :rtype: None
116    """
117# (0,0,0) -> ax
118    ax.plot([0.0, aij[0][0]], 
119            [0.0, aij[0][1]], 
120            [0.0, aij[0][2]], color = color)
121# (0,0,0) -> ay
122    ax.plot([0.0, aij[1][0]], 
123            [0.0, aij[1][1]], 
124            [0.0, aij[1][2]], color = color)
125# (0,0,0) -> az
126    ax.plot([0.0, aij[2][0]], 
127            [0.0, aij[2][1]], 
128            [0.0, aij[2][2]], color = color)
129
130# ax -> ax + ay
131    ax.plot([aij[0][0], aij[0][0] + aij[1][0]], 
132            [aij[0][1], aij[0][1] + aij[1][1]], 
133            [aij[0][2], aij[0][2] + aij[1][2]], color = color)
134# ax -> ax + az
135    ax.plot([aij[0][0], aij[0][0] + aij[2][0]], 
136            [aij[0][1], aij[0][1] + aij[2][1]], 
137            [aij[0][2], aij[0][2] + aij[2][2]], color = color)
138
139# ay -> ay + ax
140    ax.plot([aij[1][0], aij[1][0] + aij[0][0]], 
141            [aij[1][1], aij[1][1] + aij[0][1]], 
142            [aij[1][2], aij[1][2] + aij[0][2]], color = color)
143# ay -> ay + az
144    ax.plot([aij[1][0], aij[1][0] + aij[2][0]], 
145            [aij[1][1], aij[1][1] + aij[2][1]], 
146            [aij[1][2], aij[1][2] + aij[2][2]], color = color)
147
148# az -> az + ax
149    ax.plot([aij[2][0], aij[2][0] + aij[0][0]], 
150            [aij[2][1], aij[2][1] + aij[0][1]], 
151            [aij[2][2], aij[2][2] + aij[0][2]], color = color)
152# az -> ax + ay
153    ax.plot([aij[2][0], aij[2][0] + aij[1][0]], 
154            [aij[2][1], aij[2][1] + aij[1][1]], 
155            [aij[2][2], aij[2][2] + aij[1][2]], color = color)
156
157# ax + ay -> ax + ay + az
158    ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]], 
159            [aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]], 
160            [aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
161
162# ax + az -> ax + ay + az
163    ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]], 
164            [aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]], 
165            [aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
166
167# ay + az -> ax + ay + az
168    ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]], 
169            [aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]], 
170            [aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
171
172def draw_unitcell(ax, sites, aij, nrange, color = 'black'):
173    """
174    概要: 結晶の単位格子とその中の原子を3Dプロットに描画します。
175    詳細説明: draw_box関数を呼び出して単位格子を描画し、その後、sitesリスト内の原子を分数座標から
176              デカルト座標に変換してプロットします。nrangeは描画する単位格子の範囲を指定しますが、
177              このスクリプトのmain関数では現在呼び出されていません。
178    引数:
179        :param ax: matplotlibの3D軸オブジェクト。
180        :type ax: matplotlib.axes._subplots.Axes3DSubplot
181        :param sites: サイト情報のリスト。各サイトは[atom_name, site_label, atomic_number, atomic_mass, charge, radius, color, position]の形式。
182        :type sites: list of list
183        :param aij: (3, 3)のndarray、格子ベクトルa, b, cを表す。
184        :type aij: numpy.ndarray
185        :param nrange: 描画する単位格子の範囲。[[xmin, xmax], [ymin, ymax], [zmin, zmax]]の形式。
186        :type nrange: list of list of int
187        :param color: 単位格子を描画する線の色。デフォルトは'black'。
188        :type color: str
189    戻り値:
190        :returns: なし
191        :rtype: None
192    """
193    draw_box(ax, aij, nrange, color)
194
195    if sites is None:
196        return
197
198    for site in sites:
199        name, label, z, M, q, r, color, pos = site
200        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
201        for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
202         for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
203          for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
204            posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
205            if    posn[0] < nrange[0][0] or nrange[0][1] < posn[0]  \
206               or posn[1] < nrange[1][0] or nrange[1][1] < posn[1]  \
207               or posn[2] < nrange[2][0] or nrange[2][1] < posn[2]:
208                  continue
209
210            x, y, z = fractional_to_cartesian(posn, aij)
211            ax.scatter([x], [y], [z], marker = 'o', c = color, s = kr *r)
212
213
214def main():
215    """
216    概要: マデルングポテンシャルの計算と結果のプロットを実行します。
217    詳細説明:
218        格子定数から格子ベクトルや逆格子ベクトルを計算し、その情報を表示します。
219        指定された範囲rmax内で、原点にあるイオンに対するマデルングポテンシャルを
220        単純総和法で計算します。
221        計算されたポテンシャルは、距離rに対するグラフとして表示されます。
222        プログラム起動時のコマンドライン引数でrmaxとnrを設定することができます。
223    戻り値:
224        :returns: なし
225        :rtype: None
226    """
227    print("")
228    print("Lattice parameters:", lattice_parameters)
229    aij = cal_lattice_vectors(lattice_parameters)
230    print("Lattice vectors:")
231    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
232    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
233    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
234    inf = cal_metrics(lattice_parameters)
235    gij = inf['gij']
236    print("Metric tensor:")
237    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
238    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
239    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
240    volume = cal_volume(aij)
241    print("Volume: {:12.4g} A^3".format(volume))
242
243    print("")
244    print("Unit cell volume: {:12.4g} A^3".format(volume))
245    Raij  = cal_reciprocal_lattice_vectors(aij)
246    Rlatt = cal_reciprocal_lattice_parameters(Raij)
247    Rinf  = cal_metrics(Rlatt)
248    Rgij  = Rinf['gij']
249    print("Reciprocal lattice parameters:", Rlatt)
250    print("Reciprocal lattice vectors:")
251    print("  Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
252    print("  Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
253    print("  Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
254    print("Reciprocal lattice metric tensor:")
255    print("  Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
256    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
257    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
258    Rvolume = cal_volume(Raij)
259    print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
260
261# Calculate the range of unit cells
262    nxmax = int(rmax / lattice_parameters[0]) + 1
263    nymax = int(rmax / lattice_parameters[1]) + 1
264    nzmax = int(rmax / lattice_parameters[2]) + 1
265    print("")
266    print("nmax:", nxmax, nymax, nzmax)
267
268# Calculate Madelung potential around the zero-th ion
269# First store differential potential to MPdiff
270    rlist  = [rmin + i * rstep for i in range(nr)]
271    MPdiff = np.zeros(nr)
272    name0, label0, z0, M0, q0, r0, color0, pos0 = sites[0]
273    Ke = e * e / 4.0 / pi / e0                  # in MKS
274    for iz in range(-nzmax, nzmax+1):
275     for iy in range(-nymax, nymax+1):
276      for ix in range(-nxmax, nxmax+1):
277        for isite1 in range(len(sites)):
278            site1 = sites[isite1]
279            name1, label1, z1, M1, q1, r1, color1, pos1 = site1
280            r  = distance(pos0, pos1 + np.array([ix, iy, iz]), gij)
281            ir = int((r - rmin) / rstep)
282            if r < rmin or ir < 0 or nr <= ir:
283                 continue
284
285            MPdiff[ir] += Ke * q1 / (r * 1.0e-10) / e   # in eV
286
287#                print("  {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
288#                    .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))
289
290    print("")
291    print("r (A)      Madelung potential (eV)")
292    MP = np.empty(nr)
293    MP[0] = MPdiff[0]
294    print("{:10.4g}   {:12.6g}".format(rlist[0], MP[0]))
295    for i in range (1, len(MPdiff)):
296        MP[i] = MP[i-1] + MPdiff[i]
297        print("{:10.4g}   {:12.6g}".format(rlist[i], MP[i]))
298    
299    fig = plt.figure(figsize = figsize)
300    ax = fig.add_subplot(111)
301
302    ax.plot(rlist, MP)
303    ax.set_xlabel('r / angstrom')
304    ax.set_ylabel('Electrostatic potential / eV')
305
306    plt.show()
307
308    
309    terminate()
310
311
312if __name__ == '__main__':
313    argv = sys.argv
314    narg = len(argv)
315    if narg >= 2:
316        rmax = float(argv[1])
317    if narg >= 3:
318        nr = int(argv[2])
319
320    main()