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

crystal_MP_Ewald.py をダウンロード

crystal_MP_Ewald.py
crystal_MP_Ewald.py
  1"""
  2エワルド法を用いて結晶のマデルングポテンシャルを計算するスクリプトです。
  3
  4概要:
  5    このスクリプトは、指定された結晶の格子パラメータとサイト情報に基づき、
  6    エワルド法を適用してマデルングポテンシャルを計算します。
  7    実空間和、逆空間和、および自己項の3つの成分を合計して、最終的なポテンシャルを求めます。
  8
  9詳細説明:
 10    プログラムは以下の手順でマデルングポテンシャルを計算します。
 11    1.  初期化された格子パラメータとサイト情報（原子の種類、電荷、位置など）を使用します。
 12    2.  コマンドライン引数からエワルドパラメータ alpha と計算精度 prec を受け取ることができます。
 13    3.  格子ベクトル、逆格子ベクトル、体積、および関連するメトリックテンソルを計算し表示します。
 14    4.  エワルドパラメータに基づき、実空間および逆空間の計算範囲（rdmax, G2max）を決定します。
 15    5.  選択された中心サイトに対する実空間和 (UC1) を計算します。この項は実空間でのクーロン相互作用を表します。
 16    6.  逆空間和 (UC2) を計算します。この項は逆格子空間でのクーロン相互作用を表し、高速フーリエ変換に似た形式です。
 17    7.  自己相互作用項 (UC3) を計算します。これは原子自身の電場による自己エネルギーを補正する項です。
 18    8.  これら3つの項を合計し、マデルングポテンシャルおよびマデルング定数をJouleとeV単位で出力します。
 19    9.  計算にかかった時間も合わせて表示されます。
 20
 21関連リンク:
 22    crystal_MP_Ewald_usage
 23    このスクリプトは tkcrystalbase.py モジュールに定義された関数を使用します。
 24"""
 25
 26import sys
 27import os
 28import time
 29from math import erf, erfc
 30from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, log10, sqrt
 31import numpy as np
 32from numpy import linalg as la
 33import matplotlib.pyplot as plt
 34
 35from tkcrystalbase import *
 36
 37
 38pi          = 3.14159265358979323846
 39pi2         = pi + pi
 40torad       = 0.01745329251944 # rad/deg";
 41todeg       = 57.29577951472   # deg/rad";
 42basee       = 2.71828183
 43
 44h           = 6.6260755e-34    # Js";
 45h_bar       = 1.05459e-34      # "Js";
 46hbar        = h_bar
 47c           = 2.99792458e8     # m/s";
 48e           = 1.60218e-19      # C";
 49me          = 9.1093897e-31    # kg";
 50mp          = 1.6726231e-27    # kg";
 51mn          = 1.67495e-27      # kg";
 52u0          = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
 53e0          = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
 54e2_4pie0    = 2.30711e-28      # Nm<sup>2</sup>";
 55a0          = 5.29177e-11      # m";
 56kB          = 1.380658e-23     # JK<sup>-1</sup>";
 57NA          = 6.0221367e23     # mol<sup>-1</sup>";
 58R           = 8.31451          # J/K/mol";
 59F           = 96485.3          # C/mol";
 60g           = 9.81             # m/s2";
 61
 62
 63# Lattice parameters (angstrom and degree)
 64#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
 65lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
 66#lattice_parameters = [ 1.0, 1.0, 1.0, 90.0, 90.0, 90.0]
 67
 68# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
 69sites = [
 70         ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 71        ,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.5, 0.5])]
 72        ,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.0, 0.5])]
 73        ,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.5, 0.5, 0.0])]
 74        ,['Cl', 'Cl1', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
 75        ,['Cl', 'Cl2', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
 76        ,['Cl', 'Cl3', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
 77        ,['Cl', 'Cl4', 17, 35.4527,  -1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
 78        ]
 79
 80# Minimum distance to judge idential site
 81rmin = 0.1
 82
 83# Ewald alpha parameter
 84ew_alpha = 0.3
 85
 86# Precision
 87prec = 1.0e-5
 88
 89
 90def usage():
 91    """
 92    概要:
 93        プログラムの正しい使用方法を標準出力に表示します。
 94    詳細説明:
 95        プログラムをコマンドラインから実行する際の引数のフォーマットと例を示します。
 96    """
 97    print("")
 98    print("Usage: python {} alpha prec".format(argv[0]))
 99    print("   ex: python {} {} {}".format(argv[0], ew_alpha, prec))
100    print("")
101
102def terminate():
103    """
104    概要:
105        使用方法を表示した後、プログラムを終了します。
106    詳細説明:
107        usage() 関数を呼び出し、その後にプログラムを強制終了します。
108    """
109    usage()
110    exit()
111    
112
113def main():
114    """
115    概要:
116        エワルド法によりマデルングポテンシャルを計算し、結果を表示します。
117    詳細説明:
118        この関数は、設定された格子パラメータとサイト情報に基づき、
119        エワルド法を用いて結晶のマデルングポテンシャルを計算し、その結果を標準出力に表示します。
120        具体的な計算手順は以下の通りです。
121
122        1.  設定された格子パラメータから格子ベクトル、メトリックテンソル、単位胞の体積を計算し、表示します。
123            これには cal_lattice_vectors と cal_metrics 関数が使用されます。
124        2.  逆格子パラメータ、逆格子ベクトル、逆格子メトリックテンソル、逆格子単位胞の体積を計算し、表示します。
125            これには cal_reciprocal_lattice_vectors と cal_reciprocal_lattice_parameters 関数が使用されます。
126        3.  コマンドライン引数またはデフォルト設定で指定されたエワルドパラメータ ew_alpha と計算精度 prec を表示します。
127        4.  これらのパラメータに基づき、実空間和の最大距離 rdmax と逆空間和の最大Gベクトル二乗値 G2max を決定し、
128            それぞれに対応する最大繰返し回数 nrmax および hgmax を推定し、表示します。
129        5.  実空間和 UC1 の計算を実行します。これは、中心サイトと周期的に配置された他のサイトとの間のクーロン相互作用を、
130            erfc 関数を用いて収束させた合計です。
131        6.  逆空間和 UC2 の計算を実行します。これは逆格子空間における電荷分布のフーリエ成分の相互作用の合計で、
132            Gベクトルが G2max を超えない範囲で計算されます。
133        7.  自己相互作用項 UC3 を計算します。これはエワルド法の導入により生じる、
134            原子自身の電場による自己エネルギーを補正する項です。
135        8.  各計算フェーズ（実空間和、逆空間和、合計）にかかった時間を表示します。
136        9.  計算された UC1, UC2, UC3 の3つの項を合計し、最終的なマデルングポテンシャル MP を算出します。
137        10. 算出されたマデルングポテンシャル MP をJoule単位とeV単位で表示します。
138        11. 選択された中心サイトの電荷 qi と格子定数 a を用いて、マデルング定数を計算し、表示します。
139    """
140    print("")
141    print("Lattice parameters:", lattice_parameters)
142    aij = cal_lattice_vectors(lattice_parameters)
143    print("Lattice vectors:")
144    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
145    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
146    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
147    inf = cal_metrics(lattice_parameters)
148    gij = inf['gij']
149    print("Metric tensor:")
150    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
151    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
152    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
153    volume = cal_volume(aij)
154    print("Volume: {:12.4g} A^3".format(volume))
155
156    print("")
157    print("Unit cell volume: {:12.4g} A^3".format(volume))
158    Raij  = cal_reciprocal_lattice_vectors(aij)
159    Rlatt = cal_reciprocal_lattice_parameters(Raij)
160    Rinf  = cal_metrics(Rlatt)
161    Rgij  = Rinf['gij']
162    print("Reciprocal lattice parameters:", Rlatt)
163    print("Reciprocal lattice vectors:")
164    print("  Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
165    print("  Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
166    print("  Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
167    print("Reciprocal lattice metric tensor:")
168    print("  Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
169    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
170    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
171    Rvolume = cal_volume(Raij)
172    print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
173
174    nsites = len(sites)
175    
176    print("")
177    print("Ewald parameters")
178    print("  alpha:", ew_alpha)
179    norder = -log10(prec)
180    print("  precision = {} = 10^-{}".format(prec, norder))
181
182    rdmax     = (2.26 + 0.26 * norder) / ew_alpha
183    erfc_rdmax = erfc(ew_alpha * rdmax)
184    print("  RDmax = {} A, where erfc(alpha*RDmax) = {}".format(rdmax, erfc_rdmax));
185
186    lsin  = np.empty(3, dtype = float)
187    nrmax = np.empty(3, dtype = int)
188    lsin[0] = sin(torad * lattice_parameters[3])
189    lsin[1] = sin(torad * lattice_parameters[4])
190    lsin[2] = sin(torad * lattice_parameters[5])
191    nrmax[0] = int(rdmax / sqrt(gij[0][0] * lsin[1] * lsin[2])) + 1
192    nrmax[1] = int(rdmax / sqrt(gij[1][1] * lsin[2] * lsin[0])) + 1
193    nrmax[2] = int(rdmax / sqrt(gij[2][2] * lsin[0] * lsin[1])) + 1
194    print("  nrmax:", *nrmax)
195
196    cal_N = int(4.0 / 3.0 * pi * rdmax**3 / volume * nsites)
197    print("  cal_N(real):", cal_N)
198
199    G2max = ew_alpha**2 / pi**2 * (-log(prec))
200    print("  G2max:", G2max)
201    print("      exp(-pi^2 * G2max^2 / alpha^2) = ", exp(-pi**2 * G2max**2 / ew_alpha**2))
202    lsin[0] = sin(torad * Rlatt[3])
203    lsin[1] = sin(torad * Rlatt[4])
204    lsin[2] = sin(torad * Rlatt[5])
205    hgmax = np.empty(3, dtype = int)
206    hgmax[0] = int(sqrt(G2max / (Rgij[0][0] * lsin[1] * lsin[2]))) + 1
207    hgmax[1] = int(sqrt(G2max / (Rgij[1][1] * lsin[0] * lsin[2]))) + 1
208    hgmax[2] = int(sqrt(G2max / (Rgij[2][2] * lsin[0] * lsin[1]))) + 1
209    print("  hgmax:", *hgmax)
210
211    cal_N = int(4.0 / 3.0 * pi * G2max**1.5 / Rvolume * nsites)
212    print("  cal_N(reciprocal):", cal_N)
213
214    namei, labeli, zi, Mi, qi, ri, colori, pos_i = sites[0]
215
216    stime1 = time.time()
217    UC1 = 0.0
218    for iz in range(-nrmax[2], nrmax[2]+1):
219     for iy in range(-nrmax[1], nrmax[1]+1):
220      for ix in range(-nrmax[0], nrmax[0]+1):
221        for j in range(nsites):
222            namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
223            rij  = distance(pos_i, pos_j + np.array([ix, iy, iz]), gij)
224            
225            if rij < rmin:
226                 continue
227
228            erfcar = erfc(ew_alpha * rij)
229            UC1 += qj * erfcar / (rij * 1.0e-10)   # in eV
230    etime1 = time.time()
231
232    origin = np.array([0.0, 0.0, 0.0])
233    UC2 = 0.0
234    Kexp = pi * pi / ew_alpha / ew_alpha
235    Krec = 1.0 / pi / (volume * 1.0e-30)
236#    for l in range(-hgmax[2], hgmax[2]+1):
237    for l in range(0, hgmax[2]+1):
238     for k in range(-hgmax[1], hgmax[1]+1):
239      for h in range(-hgmax[0], hgmax[0]+1):
240          G2 = distance2(origin, np.array([h, k, l]), Rgij)
241          if G2 == 0.0 or G2 > G2max:
242              continue
243
244          phi_i  = pi2 * (h * pos_i[0] + k * pos_i[1] + l * pos_i[2])
245          cosphi_i = cos(phi_i)
246          sinphi_i = sin(phi_i)
247
248          cossum_j = 0.0
249          sinsum_j = 0.0
250          for j in range(nsites):
251              namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
252
253              phi_j = pi2 * (h * pos_j[0] + k * pos_j[1] + l * pos_j[2])
254              cossum_j += qj * cos(phi_j)
255              sinsum_j += qj * sin(phi_j)
256
257          fcal = cosphi_i * cossum_j + sinphi_i * sinsum_j
258          if l != 0:
259            fcal *= 2.0
260          expg = exp(-Kexp * G2) / (G2 * 1.0e+20)
261          UC2 += Krec * expg * fcal
262    etime2 = time.time()
263
264    UC3 = qi * 2.0 * (ew_alpha * 1.0e10) / sqrt(pi)
265
266    MP = UC1 + UC2 - UC3
267    etime3 = time.time()
268
269# Coefficient to calculate electrostatic potential
270    Ke = e * e / 4.0 / pi / e0
271
272    print("")
273    print("Time for real space sum     : {:6}".format(etime1 - stime1))
274    print("Time for real reciprocal sum: {:6}".format(etime2 - etime1))
275    print("Total time                  : {:6}".format(etime3 - stime1))
276    
277    print("  Madelung potential: {:12.6g} J  (= {:12.6g} + {:12.6g} - {:12.6g})".format(Ke * MP, Ke * UC1, Ke * UC2, Ke * UC3))
278    print("  Madelung potential: {:12.6g} eV (= {:12.6g} + {:12.6g} - {:12.6g})".format(Ke / e * MP, Ke / e * UC1, Ke / e * UC2, Ke / e * UC3))
279# Charge is represented by q0 to define Madelung constant
280# Lattice parameter a is represented by q0 to define Madelung constant    
281    print("  Madelung constant: {:14.8g}".format(0.5 * MP / abs(qi) * (lattice_parameters[0] * 1.0e-10)))
282
283
284    terminate()
285
286
287if __name__ == '__main__':
288    argv = sys.argv
289    narg = len(argv)
290    if narg >= 2:
291        ew_alpha = float(argv[1])
292    if narg >= 3:
293        prec = float(argv[2])
294
295    main()