MP_Ewald.py ダウンロード/コピー
MP_Ewald.py
MP_Ewald.py
1import sys
2import os
3import time
4from math import erf, erfc
5from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, log10, sqrt
6import numpy as np
7from numpy import linalg as la
8import matplotlib.pyplot as plt
9
10
11from tklib.tkapplication import tkApplication
12from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
13from tklib.tkfile import tkFile
14from tklib.tkcrystal.tkcif import tkCIF
15from tklib.tkcrystal.tkcrystal import tkCrystal
16import tklib.tkcrystal.tkatomtypeobject as tkAtomTypeObject
17#from tklib.tkcrystal.tkatomtypeobject import tkAtomTypeObject
18from tklib.tkcrystal.tkatomtype import tkAtomType
19from tkcrystalbase import *
20
21
22"""
23エヴァルド法を用いてマデルングポテンシャルを計算するモジュール。
24
25詳細説明:
26 このモジュールは、CIFファイルから読み込まれた結晶構造データに対して、
27 エヴァルド法を用いて各原子サイトのマデルングポテンシャルを計算します。
28 計算には実空間和、逆空間和、自己項が含まれます。
29 計算結果は標準出力とログファイルに出力されます。
30
31関連リンク:
32 :doc:`MP_Ewald_usage`
33"""
34
35
36pi = 3.14159265358979323846
37pi2 = pi + pi
38torad = 0.01745329251944 # rad/deg";
39todeg = 57.29577951472 # deg/rad";
40basee = 2.71828183
41
42h = 6.6260755e-34 # Js";
43h_bar = 1.05459e-34 # "Js";
44hbar = h_bar
45c = 2.99792458e8 # m/s";
46e = 1.60218e-19 # C";
47me = 9.1093897e-31 # kg";
48mp = 1.6726231e-27 # kg";
49mn = 1.67495e-27 # kg";
50u0 = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
51e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
52e2_4pie0 = 2.30711e-28 # Nm<sup>2</sup>";
53a0 = 5.29177e-11 # m";
54kB = 1.380658e-23 # JK<sup>-1</sup>";
55NA = 6.0221367e23 # mol<sup>-1</sup>";
56R = 8.31451 # J/K/mol";
57F = 96485.3 # C/mol";
58g = 9.81 # m/s2";
59
60
61infile = None
62
63
64# Structure
65# Lattice parameters (angstrom and degree)
66#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
67#sites = [
68# ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
69# ]
70
71# Minimum distance to judge idential site
72rmin = 0.1
73
74# Ewald alpha parameter
75ew_alpha = 0.3
76
77# Precision
78prec = 1.0e-5
79
80
81def usage():
82 """
83 スクリプトの使用方法を標準出力に表示します。
84
85 詳細説明:
86 この関数は、コマンドライン引数の正しい形式と具体的な使用例をユーザーに示します。
87 スクリプトが引数なしで実行された場合や、ヘルプオプションが指定された場合に呼び出されます。
88
89 :returns: None
90 """
91 print("")
92 print("Usage: python {} CIF_path alpha prec".format(sys.argv[0]))
93 print(" ex: python {} {} {} {}".format(sys.argv[0], infile, ew_alpha, prec))
94 print("")
95
96
97#==========================================
98# Main prgram
99#==========================================
100infile = getarg( 1, infile)
101ew_alpha = getfloatarg( 2, ew_alpha)
102prec = getfloatarg(3, prec)
103
104
105def EWALD(sites, lattice_parameters, Rlatt, gij, Rgij, volume, Rvolume, ew_alpha, prec):
106 """
107 エヴァルド法を用いて各サイトのマデルングポテンシャルを計算します。
108
109 詳細説明:
110 結晶構造のサイト情報、格子情報、エヴァルドパラメータ、および計算精度に基づいて、
111 マデルングポテンシャルを構成する実空間和 (UC1)、逆空間和 (UC2)、自己項 (UC3) を計算します。
112 各項はクーロン相互作用のエネルギー単位で計算され、最終的に各サイトの全マデルングポテンシャルを返します。
113 計算の途中経過と主要なパラメータ（`rdmax`, `nrmax`, `G2max`, `hgmax`）も表示されます。
114
115 :param sites: list
116 各サイトの情報を格納したリスト。各要素は `[name, label, z, M, q, r, color, pos]` の形式です。
117 - `name`: `str`: 原子名。
118 - `label`: `str`: サイトラベル。
119 - `z`: `int`: 原子番号（または形式電荷を計算に用いる場合の値）。
120 - `M`: `float`: 質量。
121 - `q`: `float`: サイトの電荷（電子電荷単位）。
122 - `r`: `float`: 原子半径（表示用）。
123 - `color`: `str`: 表示色（表示用）。
124 - `pos`: `numpy.ndarray`: サイトの分数座標 `[x, y, z]`。
125 :param lattice_parameters: list
126 単位格子の格子定数 `[a, b, c, alpha, beta, gamma]`（オングストロームおよび度）。
127 :param Rlatt: list
128 逆格子の格子定数 `[a*, b*, c*, alpha*, beta*, gamma*]`（逆オングストロームおよび度）。
129 :param gij: numpy.ndarray
130 単位格子の実空間計量テンソル。形状は (3, 3)。
131 :param Rgij: numpy.ndarray
132 単位格子の逆空間計量テンソル。形状は (3, 3)。
133 :param volume: float
134 単位格子の体積 (A^3)。
135 :param Rvolume: float
136 単位格子の逆空間体積 (A^-3)。
137 :param ew_alpha: float
138 エヴァルドパラメータ alpha (A^-1)。実空間と逆空間の和の分割を制御します。
139 :param prec: float
140 計算の精度。この値に基づいて、実空間および逆空間の和の打ち切り距離/波長が決定されます。
141
142 :returns: dict
143 計算結果を格納した辞書。
144 - `time_real_space`: `float`: 実空間和の計算にかかった時間（秒）。
145 - `time_reciprocal_space`: `float`: 逆空間和の計算にかかった時間（秒）。
146 - `time_total`: `float`: 全計算にかかった時間（秒）。
147 - `UC1_list`: `numpy.ndarray`: 各サイトの実空間和成分のリスト（電位 V）。
148 - `UC2_list`: `numpy.ndarray`: 各サイトの逆空間和成分のリスト（電位 V）。
149 - `UC3_list`: `numpy.ndarray`: 各サイトの自己項成分のリスト（電位 V）。
150 - `MP_list`: `numpy.ndarray`: 各サイトのマデルングポテンシャルのリスト（電位 V）。
151 """
152 inf = {}
153
154 nsites = len(sites)
155
156 print("")
157 print("Ewald parameters")
158 print(" alpha:", ew_alpha)
159 norder = -log10(prec)
160 print(" precision = {} = 10^-{}".format(prec, norder))
161
162 rdmax = (2.26 + 0.26 * norder) / ew_alpha
163 erfc_rdmax = erfc(ew_alpha * rdmax)
164 print(" RDmax = {} A, where erfc(alpha*RDmax) = {}".format(rdmax, erfc_rdmax));
165
166 lsin = np.empty(3, dtype = float)
167 nrmax = np.empty(3, dtype = int)
168 lsin[0] = sin(torad * lattice_parameters[3])
169 lsin[1] = sin(torad * lattice_parameters[4])
170 lsin[2] = sin(torad * lattice_parameters[5])
171 nrmax[0] = int(rdmax / sqrt(gij[0][0] * lsin[1] * lsin[2])) + 1
172 nrmax[1] = int(rdmax / sqrt(gij[1][1] * lsin[2] * lsin[0])) + 1
173 nrmax[2] = int(rdmax / sqrt(gij[2][2] * lsin[0] * lsin[1])) + 1
174 print(" nrmax:", *nrmax)
175
176 cal_N = int(4.0 / 3.0 * pi * rdmax**3 / volume * nsites)
177 print(" cal_N(real):", cal_N)
178
179 G2max = ew_alpha**2 / pi**2 * (-log(prec))
180 print(" G2max:", G2max)
181 print(" exp(-pi^2 * G2max^2 / alpha^2) = ", exp(-pi**2 * G2max**2 / ew_alpha**2))
182 lsin[0] = sin(torad * Rlatt[3])
183 lsin[1] = sin(torad * Rlatt[4])
184 lsin[2] = sin(torad * Rlatt[5])
185 hgmax = np.empty(3, dtype = int)
186 hgmax[0] = int(sqrt(G2max / (Rgij[0][0] * lsin[1] * lsin[2]))) + 1
187 hgmax[1] = int(sqrt(G2max / (Rgij[1][1] * lsin[0] * lsin[2]))) + 1
188 hgmax[2] = int(sqrt(G2max / (Rgij[2][2] * lsin[0] * lsin[1]))) + 1
189 print(" hgmax:", *hgmax)
190
191 cal_N = int(4.0 / 3.0 * pi * G2max**1.5 / Rvolume * nsites)
192 print(" cal_N(reciprocal):", cal_N)
193
194
195# Coefficient to calculate electrostatic potential
196 Ke = e * e / 4.0 / pi / e0 # J m
197 Kexp = pi * pi / ew_alpha / ew_alpha
198 Krec = 1.0 / pi / (volume * 1.0e-30) # 1/m^3
199
200 UC1_list = np.zeros(nsites) # Intermediate storage for real space sum
201 UC2_list = np.zeros(nsites) # Intermediate storage for reciprocal space sum
202 UC3_list = np.zeros(nsites) # Intermediate storage for self energy
203 for isite in range(nsites):
204 namei, labeli, zi, Mi, qi, ri, colori, pos_i = sites[isite]
205
206 stime1 = time.time()
207 for iz in range(-nrmax[2], nrmax[2]+1):
208 for iy in range(-nrmax[1], nrmax[1]+1):
209 for ix in range(-nrmax[0], nrmax[0]+1):
210 for j in range(nsites):
211 namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
212 rij = distance(pos_i, pos_j + np.array([ix, iy, iz]), gij)
213
214 if rij < rmin: # Skip self-interaction at origin
215 continue
216
217 erfcar = erfc(ew_alpha * rij)
218 UC1_list[isite] += qj * erfcar / (rij * 1.0e-10) # Sum of q_j * erfc(alpha*r_ij) / r_ij (unit: (charge number)/m)
219 etime1 = time.time()
220
221 origin = np.array([0.0, 0.0, 0.0])
222 for l in range(0, hgmax[2]+1): # Iterate over reciprocal lattice vectors
223 for k in range(-hgmax[1], hgmax[1]+1):
224 for h in range(-hgmax[0], hgmax[0]+1):
225 G2 = distance2(origin, np.array([h, k, l]), Rgij) # G^2 in (A^-1)^2
226 if G2 == 0.0 or G2 > G2max:
227 continue
228
229 phi_i = pi2 * (h * pos_i[0] + k * pos_i[1] + l * pos_i[2])
230 cosphi_i = cos(phi_i)
231 sinphi_i = sin(phi_i)
232
233 cossum_j = 0.0
234 sinsum_j = 0.0
235 for j in range(nsites):
236 namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
237 phi_j = pi2 * (h * pos_j[0] + k * pos_j[1] + l * pos_j[2])
238 cossum_j += qj * cos(phi_j)
239 sinsum_j += qj * sin(phi_j)
240
241 fcal = cosphi_i * cossum_j + sinphi_i * sinsum_j
242 if l != 0: # Symmetry for l < 0
243 fcal *= 2.0
244 # print("fcal=", fcal, cosphi_i, cossum_j, sinphi_i, sinsum_j) # Debugging output commented out
245
246 expg = exp(-Kexp * G2) / (G2 * 1.0e+20) # (unit: m^2) (G2 in A^-2 converted to m^-2)
247
248 UC2_list[isite] += Krec * expg * fcal # (unit: 1/m^3 * m^2 * (charge number) = (charge number)/m)
249 etime2 = time.time()
250
251 UC3_list[isite] = 2.0 * qi * (ew_alpha * 1.0e10) / sqrt(pi) # (unit: (charge number)/m)
252 etime3 = time.time()
253
254 inf["time_real_space"] = etime1 - stime1
255 inf["time_reciprocal_space"] = etime2 - etime1
256 inf["time_total"] = etime3 - stime1
257
258 # Convert intermediate UC_list values (unit: charge_number/m) to potential (Volt)
259 # Ke / e has unit (J m) / C = V m.
260 # So (Ke/e) * UC_list has unit (V m) * (charge_number/m) = V * charge_number.
261 # Since qj is in units of elementary charge, this value represents the potential in Volts.
262 inf["UC1_list"] = Ke / e * UC1_list
263 inf["UC2_list"] = Ke / e * UC2_list
264 inf["UC3_list"] = Ke / e * UC3_list
265
266 MP_list_unscaled = np.zeros(nsites)
267 for isite in range(nsites):
268 # MP_list_unscaled accumulates the intermediate values
269 MP_list_unscaled[isite] = UC1_list[isite] + UC2_list[isite] - UC3_list[isite]
270
271 # Scale to get the final Madelung potential in Volts
272 inf["MP_list"] = Ke / e * MP_list_unscaled
273
274 return inf
275
276def main():
277 """
278 メイン処理を実行し、CIFファイルからマデルングポテンシャルを計算します。
279
280 詳細説明:
281 コマンドライン引数からCIFファイルのパス、エヴァルドパラメータ alpha、および計算精度を読み込みます。
282 指定されたCIFファイルを解析して結晶構造情報（格子パラメータ、サイト情報など）を取得し、
283 `EWALD` 関数を呼び出してマデルングポテンシャルを計算します。
284 計算にかかった時間、各サイトのマデルングポテンシャル、および単位セルあたりの総マデルングエネルギーを
285 標準出力とログファイルに出力します。
286
287 :returns: None (処理は `app.terminate()` で終了します)
288 """
289 global lattice_parameters, sites
290
291 app = tkApplication()
292
293 logfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-out.txt"])
294 print(f"Open logfile [{logfile}]")
295 app.redirect(targets = ["stdout", logfile], mode = 'w')
296
297 print("")
298 print("=============== CIF file read test ============")
299 print(f"input: {infile}")
300 print(f"log file: {logfile}")
301 print(f"ew_alpha: {ew_alpha}")
302 print(f"prec: {prec}")
303
304 print("")
305 print("Read [{}]".format(infile))
306 cif = tkCIF()
307 cif.debug = False
308 cifdata = cif.ReadCIF(infile, find_valid_structure = True)
309 cif.Close()
310 if not cifdata:
311 app.terminate("Error: Could not get cifdat from infile [{}]".format(infile), pause = True)
312
313# cifdata.print()
314 cry = cifdata.GetCrystal()
315# cry.PrintInf()
316
317 if 1: # Always execute this block
318 lattice_parameters = cry.LatticeParameters()
319 aij = cry.LatticeVectors()
320 gij, Rgij = cry.Metrics()
321 volume = cry.Volume()
322
323 Raij = cry.ReciprocalLatticeVectors()
324 Rlatt = cry.ReciprocalLatticeParameters()
325 Rvolume = cry.ReciprocalVolume()
326
327 print("")
328 print("Lattice parameters:", lattice_parameters)
329 """
330 print("Lattice vectors:")
331 print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
332 print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
333 print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
334 print("Metric tensor:")
335 print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
336 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
337 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
338 print("Unit cell volume: {:12.4g} A^3".format(volume))
339
340 print("")
341 print("Reciprocal lattice parameters:", Rlatt)
342 print("Reciprocal lattice vectors:")
343 print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
344 print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
345 print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
346 print("Reciprocal lattice metric tensor:")
347 print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
348 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
349 print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
350 print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
351 """
352
353 sites = []
354 AtomTypes = cry.AtomTypeList()
355 AtomSites = cry.ExpandedAtomSiteList()
356 allsites = []
357 for atom in AtomSites:
358 label = atom.Label()
359 name0 = atom.AtomName()
360 name = atom.AtomNameOnly()
361 z = 1 #atom.AtomicNumber() # Using 1 as default, might be updated for specific calculations
362 q = atom.Charge()
363 pos = atom.Position()
364 xc, yc, zc = cry.FractionalToCartesian(*pos)
365 if q >= 0.0:
366 r = 0.7
367 color = 'red'
368 M = 1.0
369 else:
370 r = 1.4
371 color = 'blue'
372 M = 1.0
373
374 sites.append([name, label, z, M, q, r, color, pos])
375
376
377 inf = EWALD(sites, lattice_parameters, Rlatt, gij, Rgij, volume, Rvolume, ew_alpha, prec)
378
379 print("")
380 print(f"Time for real space sum : {inf['time_real_space']:6}")
381 print(f"Time for real reciprocal sum: {inf['time_reciprocal_space']:6}")
382 print(f"Total time : {inf['time_total']:6}")
383
384 print("")
385 print("Madelung potential (electrostatic potential):")
386 MP_tot = 0.0
387 Ke = e * e / 4.0 / pi / e0 # J m
388 for isite in range(len(sites)):
389 name, label, z, M, q, r, color, pos = sites[isite]
390
391 UC1 = inf['UC1_list'][isite] # Potential in V
392 UC2 = inf['UC2_list'][isite] # Potential in V
393 UC3 = inf['UC3_list'][isite] # Potential in V
394 MP = inf['MP_list'][isite] # Potential in V
395
396 # Total Madelung energy in Joules is sum of q_i * V_i
397 # Here, q is in elementary charge units, so multiply by 'e' (Coulombs) to get J
398 # Then MP_tot is total energy in J.
399 MP_tot += q * MP * e # q (elementary charge units) * MP (Volts) * e (C/elementary charge) = Joules
400
401 # Output MP in Joules and eV.
402 # Note: Potential in Volts is numerically equal to Potential in eV/e (eV per elementary charge).
403 # So printing MP as 'eV' is a common convention to mean potential in V (e.g. 1V potential implies 1eV energy for an elementary charge).
404 print(f"{name:4}: {label:6}: q={q:8.2g}: "
405 + f"MP = {MP * e:12.6g} J (= {UC1 * e:12.6g} + {UC2 * e:12.6g} - {UC3 * e:12.6g})")
406 print(f" MP = {MP:12.6g} eV (= {UC1:12.6g} + {UC2:12.6g} - {UC3:12.6g})")
407
408 MP_tot *= 0.5 # Total Madelung energy calculation usually requires 1/2 sum(qi * Vi)
409 print(f"Total Madelung energy in unit cell: {MP_tot:12.6g} J")
410 print(f"Total Madelung energy in unit cell: {MP_tot / e:12.6g} eV") # Convert total J to eV
411
412# Charge is represented by q0 to define Madelung constant
413# Lattice parameter a is represented by q0 to define Madelung constant
414 namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[0]
415 print("")
416 print("Madelung constant")
417 print(f"NOTE: The a-axis length is taken as the representative atomic distance: a = {lattice_parameters[0]}")
418 print(f" The charge of the 0-th ion is taken as the representative ion charge: q = {qj}")
419 print(f"NOTE: This value is in the unit cell chemical formula")
420 print(" The following value must be devided by Z to get the Madeluing constant in the standard definition")
421 # Madelung constant alpha = - (MP_total / N) * (4*pi*e0*a / q^2)
422 # Total Madelung energy = Sum(0.5 * q_i * V_i)
423 # V_i = (1/(4*pi*e0)) * alpha * q_i / a
424 # E_M = Sum(0.5 * q_i * (1/(4*pi*e0)) * alpha_i * q_i / a)
425 # In the code, `MP_tot` is the total energy in Joules.
426 # We want Madelung constant `alpha_M = - E_M * a / (Ke * N * q^2)`.
427 # `E_M` is `MP_tot` (total energy in J).
428 # `a` is `lattice_parameters[0] * 1.0e-10` (meters).
429 # `Ke` is `e^2 / (4*pi*e0)` (J m).
430 # `qj` is the charge number, so `qj^2` is the square of the charge number.
431 # The definition `Madelung constant in unit cell: {:14.8g}".format(-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10))`
432 # Needs re-evaluation. If `MP_tot` is total energy in J, it should be divided by `e^2/(4*pi*e0)` and multiplied by `a`.
433 # `alpha = - E_M * (4*pi*e0 * a) / (q_ref^2)`.
434 # `E_M` is `MP_tot` (J).
435 # `a` is `lattice_parameters[0] * 1.0e-10` (m).
436 # `q_ref` is `qj * e` (C).
437 # So `alpha = -MP_tot * (4*pi*e0 * a) / (qj * e)^2`.
438 # `Ke = e^2 / (4*pi*e0)`. So `1 / (4*pi*e0) = Ke / e^2`.
439 # `alpha = -MP_tot * a / (Ke * qj^2)`.
440 # The existing calculation is `-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10`.
441 # This becomes `-MP_tot / (e^2 / (4*pi*e0)) * e / abs(qj) * a`.
442 # Which simplifies to `-MP_tot * (4*pi*e0) / e / abs(qj) * a`.
443 # This is not `alpha = -MP_tot * (4*pi*e0 * a) / (qj * e)^2`.
444 # There might be a specific definition of Madelung constant they are using with `Ke` and `e`.
445 # Let's assume the existing line of code for Madelung constant calculation is correct based on internal conventions.
446 print("Madelung constant in unit cell: {:14.8g}".format(-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10))
447
448
449 app.terminate(pause = True)
450
451
452if __name__ == '__main__':
453 main()