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

crystal_MP_Evjen.py をダウンロード

crystal_MP_Evjen.py
crystal_MP_Evjen.py
  1"""
  2概要: Evjen法を用いて結晶のマデルングポテンシャルを計算します。
  3詳細説明:
  4    このスクリプトは、指定された結晶構造（格子定数とサイト情報）に基づいて、Evjen法を用いてマデルングポテンシャルおよびマデルング定数を算出します。結晶学的な各種パラメータ（格子ベクトル、体積、逆格子など）も計算・表示します。
  5    マデルングポテンシャルの計算では、中心イオンの周りに配置された他のイオンからのクーロン相互作用を、一定の範囲内でセルを拡大しながら合計します。このスクリプトはtkcrystalbase.pyモジュールに依存しています。
  6関連リンク:
  7    tkcrystalbase.py (結晶構造の基本機能を提供するモジュール)
  8"""
  9import sys
 10import os
 11from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
 12import numpy as np
 13from numpy import linalg as la
 14import matplotlib.pyplot as plt
 15
 16from tkcrystalbase import *
 17
 18
 19pi          = 3.14159265358979323846
 20pi2         = pi + pi
 21torad       = 0.01745329251944 # rad/deg";
 22todeg       = 57.29577951472   # deg/rad";
 23basee       = 2.71828183
 24
 25h           = 6.6260755e-34    # Js";
 26h_bar       = 1.05459e-34      # "Js";
 27hbar        = h_bar
 28c           = 2.99792458e8     # m/s";
 29e           = 1.60218e-19      # C";
 30me          = 9.1093897e-31    # kg";
 31mp          = 1.6726231e-27    # kg";
 32mn          = 1.67495e-27      # kg";
 33u0          = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
 34e0          = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
 35e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
 36a0          = 5.29177e-11      # m";
 37kB          = 1.380658e-23     # JK<sup>-1</sup>";
 38NA          = 6.0221367e23     # mol<sup>-1</sup>";
 39R           = 8.31451          # J/K/mol";
 40F           = 96485.3          # C/mol";
 41g           = 9.81             # m/s2";
 42
 43
 44# Lattice parameters (angstrom and degree)
 45#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
 46lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
 47#lattice_parameters = [ 1.0, 1.0, 1.0, 90.0, 90.0, 90.0]
 48
 49# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
 50sites = [
 51         ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 52        ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.5, 0.5])]
 53        ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.0, 0.5])]
 54        ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.0])]
 55        ,['Cl', 'Cl1', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
 56        ,['Cl', 'Cl2', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
 57        ,['Cl', 'Cl3', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
 58        ,['Cl', 'Cl4', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
 59        ]
 60
 61# Range of unit cells to regard as in the unit cell
 62nrange0 = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]]
 63
 64# Range of unit cells to calculate Madelung potential
 65nmax = 1
 66
 67# Minimum distance to judge idential site
 68rmin = 0.1
 69
 70
 71def usage():
 72    """
 73    概要: スクリプトの正しい使用方法を標準出力に表示します。
 74    詳細説明: この関数は、プログラムの引数が不適切である場合などに呼び出され、ユーザーに引数の指定方法を案内します。
 75    引数:
 76        なし
 77    戻り値:
 78        なし
 79    """
 80    print("")
 81    print("Usage: python {} rmax".format(argv[0]))
 82    print("   ex: python {} {}".format(argv[0], nmax))
 83    print("")
 84
 85def terminate():
 86    """
 87    概要: usage関数を呼び出して使用方法を表示し、プログラムを終了します。
 88    詳細説明: 通常、コマンドライン引数のエラー時など、プログラムの実行を中断する必要がある場合に呼び出されます。
 89    引数:
 90        なし
 91    戻り値:
 92        なし
 93    """
 94    usage()
 95    exit()
 96
 97
 98def main():
 99    """
100    概要: プログラムの主要な処理を実行します。
101    詳細説明:
102        格子定数から格子ベクトル、体積、逆格子情報を計算し、表示します。
103        次に、単位セル内のサイトの多重度をカウントし、全てのサイト情報を拡張します。
104        最後に、Evjen法を用いて指定されたnmaxの範囲で中心イオン周りのマデルングポテンシャルを計算し、その結果からマデルングポテンシャル（J、eV単位）とマデルング定数を計算・表示します。
105        nmaxはコマンドライン引数から取得されます。
106    引数:
107        なし
108    戻り値:
109        なし
110    """
111    print("")
112    print("Lattice parameters:", lattice_parameters)
113    aij = cal_lattice_vectors(lattice_parameters)
114    print("Lattice vectors:")
115    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
116    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
117    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
118    inf = cal_metrics(lattice_parameters)
119    gij = inf['gij']
120    print("Metric tensor:")
121    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
122    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
123    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
124    volume = cal_volume(aij)
125    print("Volume: {:12.4g} A^3".format(volume))
126
127    print("")
128    print("Unit cell volume: {:12.4g} A^3".format(volume))
129    Raij  = cal_reciprocal_lattice_vectors(aij)
130    Rlatt = cal_reciprocal_lattice_parameters(Raij)
131    Rinf  = cal_metrics(Rlatt)
132    Rgij  = Rinf['gij']
133    print("Reciprocal lattice parameters:", Rlatt)
134    print("Reciprocal lattice vectors:")
135    print("  Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
136    print("  Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
137    print("  Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
138    print("Reciprocal lattice metric tensor:")
139    print("  Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
140    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
141    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
142    Rvolume = cal_volume(Raij)
143    print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
144
145# Count multiplicity in the unit cell
146    mult = np.zeros(len(sites), dtype = int)
147    extpos = []
148    for isite in range(len(sites)):
149        name, label, z, M, q, r, color, pos = sites[isite]
150        pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
151        for iz in range(int(nrange0[2][0]) - 1, int(nrange0[2][1]) + 1):
152         for iy in range(int(nrange0[1][0]) - 1, int(nrange0[1][1]) + 1):
153          for ix in range(int(nrange0[0][0]) - 1, int(nrange0[0][1]) + 1):
154            posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
155            if    posn[0] < nrange0[0][0] or nrange0[0][1] < posn[0]  \
156               or posn[1] < nrange0[1][0] or nrange0[1][1] < posn[1]  \
157               or posn[2] < nrange0[2][0] or nrange0[2][1] < posn[2]:
158                  continue
159
160            mult[isite] += 1
161            extpos.append([name, label, z, M, q, r, [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz], isite])
162
163    print("")
164    print("Site information (all sites in the unit cell with the range:", nrange0, ")")
165    qtot = 0.0
166    for isite in range(len(extpos)):
167        name, label, z, M, q, r, pos, isite0 = extpos[isite]
168        m = mult[isite0]
169        w = 1.0 / m
170        print("  {:4} ({:8.3g}, {:8.3g}, {:8.3g}) q={:5.3g} mult={:2d} weight={:8.4g}".format(label, *pos, q, m, w))
171        qtot += q * w
172    print("qtot=", qtot)
173
174    print("")
175    print("Calculate Madelung potential around the zero-th ion by Evjen method")
176    print("  nmax:", nmax)
177    name0, label0, z0, M0, q0, r0, pos0, isite = extpos[0]
178    print("  Origin: {} ({}, {}, {})".format(label0, *pos0))
179    MP = 0.0
180    for iz in range(-nmax, nmax):
181     for iy in range(-nmax, nmax):
182      for ix in range(-nmax, nmax):
183        for isite1 in range(len(extpos)):
184            extsite1 = extpos[isite1]
185            name1, label1, z1, M1, q1, r1, pos1, idxsite = extsite1
186            w = 1.0 / mult[idxsite]
187
188            r = distance(pos0, pos1 + np.array([ix, iy, iz]), gij)
189            if r < rmin:
190                 continue
191
192            p = q1 / (r * 1.0e-10) * w
193            MP += p
194#            print(" ({:2d},{:2d},{:2d})+({:8.3g}, {:8.3g}, {:8.3g}) d={:8.3g} q={:2g} w={:8.4g} p={:8.4g} MP={:8.4g}".format(ix, iy, iz, *pos1, r, q1, w, p, MP))
195
196# Coefficient to calculate electrostatic potential
197    Ke = e * e / 4.0 / pi / e0
198
199    print("")
200    print("  Madelung potential: {:12.6g} J".format(Ke * MP))
201    print("  Madelung potential: {:12.6g} eV".format(Ke / e * MP))
202# Charge is represented by q0 to define Madelung constant
203# Lattice parameter a is represented by q0 to define Madelung constant    
204    print("  Madelung constant: {:14.8g}".format(0.5 * MP / abs(q0) * (lattice_parameters[0] * 1.0e-10)))
205
206
207    terminate()
208
209
210if __name__ == '__main__':
211    argv = sys.argv
212    narg = len(argv)
213    if narg >= 2:
214        nmax = int(argv[1])
215
216    main()