import sys
import os
import time
from math import erf, erfc
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, log10, sqrt
import numpy as np
from numpy import linalg as la
import matplotlib.pyplot as plt
from tklib.tkapplication import tkApplication
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tkfile import tkFile
from tklib.tkcrystal.tkcif import tkCIF
from tklib.tkcrystal.tkcrystal import tkCrystal
import tklib.tkcrystal.tkatomtypeobject as tkAtomTypeObject
#from tklib.tkcrystal.tkatomtypeobject import tkAtomTypeObject
from tklib.tkcrystal.tkatomtype import tkAtomType
from tkcrystalbase import *
"""
エヴァルド法を用いてマデルングポテンシャルを計算するモジュール。
詳細説明:
このモジュールは、CIFファイルから読み込まれた結晶構造データに対して、
エヴァルド法を用いて各原子サイトのマデルングポテンシャルを計算します。
計算には実空間和、逆空間和、自己項が含まれます。
計算結果は標準出力とログファイルに出力されます。
関連リンク:
:doc:`MP_Ewald_usage`
"""
pi = 3.14159265358979323846
pi2 = pi + pi
torad = 0.01745329251944 # rad/deg";
todeg = 57.29577951472 # deg/rad";
basee = 2.71828183
h = 6.6260755e-34 # Js";
h_bar = 1.05459e-34 # "Js";
hbar = h_bar
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
me = 9.1093897e-31 # kg";
mp = 1.6726231e-27 # kg";
mn = 1.67495e-27 # kg";
u0 = 4.0 * 3.14*1e-7; # . "Ns<sup>2</sup>C<sup>-2</sup>";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e2_4pie0 = 2.30711e-28 # Nm<sup>2</sup>";
a0 = 5.29177e-11 # m";
kB = 1.380658e-23 # JK<sup>-1</sup>";
NA = 6.0221367e23 # mol<sup>-1</sup>";
R = 8.31451 # J/K/mol";
F = 96485.3 # C/mol";
g = 9.81 # m/s2";
infile = None
# Structure
# Lattice parameters (angstrom and degree)
#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
#sites = [
# ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
# ]
# Minimum distance to judge idential site
rmin = 0.1
# Ewald alpha parameter
ew_alpha = 0.3
# Precision
prec = 1.0e-5
[ドキュメント]
def usage():
"""
スクリプトの使用方法を標準出力に表示します。
詳細説明:
この関数は、コマンドライン引数の正しい形式と具体的な使用例をユーザーに示します。
スクリプトが引数なしで実行された場合や、ヘルプオプションが指定された場合に呼び出されます。
:returns: None
"""
print("")
print("Usage: python {} CIF_path alpha prec".format(sys.argv[0]))
print(" ex: python {} {} {} {}".format(sys.argv[0], infile, ew_alpha, prec))
print("")
#==========================================
# Main prgram
#==========================================
infile = getarg( 1, infile)
ew_alpha = getfloatarg( 2, ew_alpha)
prec = getfloatarg(3, prec)
[ドキュメント]
def EWALD(sites, lattice_parameters, Rlatt, gij, Rgij, volume, Rvolume, ew_alpha, prec):
"""
エヴァルド法を用いて各サイトのマデルングポテンシャルを計算します。
詳細説明:
結晶構造のサイト情報、格子情報、エヴァルドパラメータ、および計算精度に基づいて、
マデルングポテンシャルを構成する実空間和 (UC1)、逆空間和 (UC2)、自己項 (UC3) を計算します。
各項はクーロン相互作用のエネルギー単位で計算され、最終的に各サイトの全マデルングポテンシャルを返します。
計算の途中経過と主要なパラメータ（`rdmax`, `nrmax`, `G2max`, `hgmax`）も表示されます。
:param sites: list
各サイトの情報を格納したリスト。各要素は `[name, label, z, M, q, r, color, pos]` の形式です。
- `name`: `str`: 原子名。
- `label`: `str`: サイトラベル。
- `z`: `int`: 原子番号（または形式電荷を計算に用いる場合の値）。
- `M`: `float`: 質量。
- `q`: `float`: サイトの電荷（電子電荷単位）。
- `r`: `float`: 原子半径（表示用）。
- `color`: `str`: 表示色（表示用）。
- `pos`: `numpy.ndarray`: サイトの分数座標 `[x, y, z]`。
:param lattice_parameters: list
単位格子の格子定数 `[a, b, c, alpha, beta, gamma]`（オングストロームおよび度）。
:param Rlatt: list
逆格子の格子定数 `[a*, b*, c*, alpha*, beta*, gamma*]`（逆オングストロームおよび度）。
:param gij: numpy.ndarray
単位格子の実空間計量テンソル。形状は (3, 3)。
:param Rgij: numpy.ndarray
単位格子の逆空間計量テンソル。形状は (3, 3)。
:param volume: float
単位格子の体積 (A^3)。
:param Rvolume: float
単位格子の逆空間体積 (A^-3)。
:param ew_alpha: float
エヴァルドパラメータ alpha (A^-1)。実空間と逆空間の和の分割を制御します。
:param prec: float
計算の精度。この値に基づいて、実空間および逆空間の和の打ち切り距離/波長が決定されます。
:returns: dict
計算結果を格納した辞書。
- `time_real_space`: `float`: 実空間和の計算にかかった時間（秒）。
- `time_reciprocal_space`: `float`: 逆空間和の計算にかかった時間（秒）。
- `time_total`: `float`: 全計算にかかった時間（秒）。
- `UC1_list`: `numpy.ndarray`: 各サイトの実空間和成分のリスト（電位 V）。
- `UC2_list`: `numpy.ndarray`: 各サイトの逆空間和成分のリスト（電位 V）。
- `UC3_list`: `numpy.ndarray`: 各サイトの自己項成分のリスト（電位 V）。
- `MP_list`: `numpy.ndarray`: 各サイトのマデルングポテンシャルのリスト（電位 V）。
"""
inf = {}
nsites = len(sites)
print("")
print("Ewald parameters")
print(" alpha:", ew_alpha)
norder = -log10(prec)
print(" precision = {} = 10^-{}".format(prec, norder))
rdmax = (2.26 + 0.26 * norder) / ew_alpha
erfc_rdmax = erfc(ew_alpha * rdmax)
print(" RDmax = {} A, where erfc(alpha*RDmax) = {}".format(rdmax, erfc_rdmax));
lsin = np.empty(3, dtype = float)
nrmax = np.empty(3, dtype = int)
lsin[0] = sin(torad * lattice_parameters[3])
lsin[1] = sin(torad * lattice_parameters[4])
lsin[2] = sin(torad * lattice_parameters[5])
nrmax[0] = int(rdmax / sqrt(gij[0][0] * lsin[1] * lsin[2])) + 1
nrmax[1] = int(rdmax / sqrt(gij[1][1] * lsin[2] * lsin[0])) + 1
nrmax[2] = int(rdmax / sqrt(gij[2][2] * lsin[0] * lsin[1])) + 1
print(" nrmax:", *nrmax)
cal_N = int(4.0 / 3.0 * pi * rdmax**3 / volume * nsites)
print(" cal_N(real):", cal_N)
G2max = ew_alpha**2 / pi**2 * (-log(prec))
print(" G2max:", G2max)
print(" exp(-pi^2 * G2max^2 / alpha^2) = ", exp(-pi**2 * G2max**2 / ew_alpha**2))
lsin[0] = sin(torad * Rlatt[3])
lsin[1] = sin(torad * Rlatt[4])
lsin[2] = sin(torad * Rlatt[5])
hgmax = np.empty(3, dtype = int)
hgmax[0] = int(sqrt(G2max / (Rgij[0][0] * lsin[1] * lsin[2]))) + 1
hgmax[1] = int(sqrt(G2max / (Rgij[1][1] * lsin[0] * lsin[2]))) + 1
hgmax[2] = int(sqrt(G2max / (Rgij[2][2] * lsin[0] * lsin[1]))) + 1
print(" hgmax:", *hgmax)
cal_N = int(4.0 / 3.0 * pi * G2max**1.5 / Rvolume * nsites)
print(" cal_N(reciprocal):", cal_N)
# Coefficient to calculate electrostatic potential
Ke = e * e / 4.0 / pi / e0 # J m
Kexp = pi * pi / ew_alpha / ew_alpha
Krec = 1.0 / pi / (volume * 1.0e-30) # 1/m^3
UC1_list = np.zeros(nsites) # Intermediate storage for real space sum
UC2_list = np.zeros(nsites) # Intermediate storage for reciprocal space sum
UC3_list = np.zeros(nsites) # Intermediate storage for self energy
for isite in range(nsites):
namei, labeli, zi, Mi, qi, ri, colori, pos_i = sites[isite]
stime1 = time.time()
for iz in range(-nrmax[2], nrmax[2]+1):
for iy in range(-nrmax[1], nrmax[1]+1):
for ix in range(-nrmax[0], nrmax[0]+1):
for j in range(nsites):
namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
rij = distance(pos_i, pos_j + np.array([ix, iy, iz]), gij)
if rij < rmin: # Skip self-interaction at origin
continue
erfcar = erfc(ew_alpha * rij)
UC1_list[isite] += qj * erfcar / (rij * 1.0e-10) # Sum of q_j * erfc(alpha*r_ij) / r_ij (unit: (charge number)/m)
etime1 = time.time()
origin = np.array([0.0, 0.0, 0.0])
for l in range(0, hgmax[2]+1): # Iterate over reciprocal lattice vectors
for k in range(-hgmax[1], hgmax[1]+1):
for h in range(-hgmax[0], hgmax[0]+1):
G2 = distance2(origin, np.array([h, k, l]), Rgij) # G^2 in (A^-1)^2
if G2 == 0.0 or G2 > G2max:
continue
phi_i = pi2 * (h * pos_i[0] + k * pos_i[1] + l * pos_i[2])
cosphi_i = cos(phi_i)
sinphi_i = sin(phi_i)
cossum_j = 0.0
sinsum_j = 0.0
for j in range(nsites):
namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[j]
phi_j = pi2 * (h * pos_j[0] + k * pos_j[1] + l * pos_j[2])
cossum_j += qj * cos(phi_j)
sinsum_j += qj * sin(phi_j)
fcal = cosphi_i * cossum_j + sinphi_i * sinsum_j
if l != 0: # Symmetry for l < 0
fcal *= 2.0
# print("fcal=", fcal, cosphi_i, cossum_j, sinphi_i, sinsum_j) # Debugging output commented out
expg = exp(-Kexp * G2) / (G2 * 1.0e+20) # (unit: m^2) (G2 in A^-2 converted to m^-2)
UC2_list[isite] += Krec * expg * fcal # (unit: 1/m^3 * m^2 * (charge number) = (charge number)/m)
etime2 = time.time()
UC3_list[isite] = 2.0 * qi * (ew_alpha * 1.0e10) / sqrt(pi) # (unit: (charge number)/m)
etime3 = time.time()
inf["time_real_space"] = etime1 - stime1
inf["time_reciprocal_space"] = etime2 - etime1
inf["time_total"] = etime3 - stime1
# Convert intermediate UC_list values (unit: charge_number/m) to potential (Volt)
# Ke / e has unit (J m) / C = V m.
# So (Ke/e) * UC_list has unit (V m) * (charge_number/m) = V * charge_number.
# Since qj is in units of elementary charge, this value represents the potential in Volts.
inf["UC1_list"] = Ke / e * UC1_list
inf["UC2_list"] = Ke / e * UC2_list
inf["UC3_list"] = Ke / e * UC3_list
MP_list_unscaled = np.zeros(nsites)
for isite in range(nsites):
# MP_list_unscaled accumulates the intermediate values
MP_list_unscaled[isite] = UC1_list[isite] + UC2_list[isite] - UC3_list[isite]
# Scale to get the final Madelung potential in Volts
inf["MP_list"] = Ke / e * MP_list_unscaled
return inf
[ドキュメント]
def main():
"""
メイン処理を実行し、CIFファイルからマデルングポテンシャルを計算します。
詳細説明:
コマンドライン引数からCIFファイルのパス、エヴァルドパラメータ alpha、および計算精度を読み込みます。
指定されたCIFファイルを解析して結晶構造情報（格子パラメータ、サイト情報など）を取得し、
`EWALD` 関数を呼び出してマデルングポテンシャルを計算します。
計算にかかった時間、各サイトのマデルングポテンシャル、および単位セルあたりの総マデルングエネルギーを
標準出力とログファイルに出力します。
:returns: None (処理は `app.terminate()` で終了します)
"""
global lattice_parameters, sites
app = tkApplication()
logfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-out.txt"])
print(f"Open logfile [{logfile}]")
app.redirect(targets = ["stdout", logfile], mode = 'w')
print("")
print("=============== CIF file read test ============")
print(f"input: {infile}")
print(f"log file: {logfile}")
print(f"ew_alpha: {ew_alpha}")
print(f"prec: {prec}")
print("")
print("Read [{}]".format(infile))
cif = tkCIF()
cif.debug = False
cifdata = cif.ReadCIF(infile, find_valid_structure = True)
cif.Close()
if not cifdata:
app.terminate("Error: Could not get cifdat from infile [{}]".format(infile), pause = True)
# cifdata.print()
cry = cifdata.GetCrystal()
# cry.PrintInf()
if 1: # Always execute this block
lattice_parameters = cry.LatticeParameters()
aij = cry.LatticeVectors()
gij, Rgij = cry.Metrics()
volume = cry.Volume()
Raij = cry.ReciprocalLatticeVectors()
Rlatt = cry.ReciprocalLatticeParameters()
Rvolume = cry.ReciprocalVolume()
print("")
print("Lattice parameters:", lattice_parameters)
"""
print("Lattice vectors:")
print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
print("Metric tensor:")
print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
print("Unit cell volume: {:12.4g} A^3".format(volume))
print("")
print("Reciprocal lattice parameters:", Rlatt)
print("Reciprocal lattice vectors:")
print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
print("Reciprocal lattice metric tensor:")
print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
"""
sites = []
AtomTypes = cry.AtomTypeList()
AtomSites = cry.ExpandedAtomSiteList()
allsites = []
for atom in AtomSites:
label = atom.Label()
name0 = atom.AtomName()
name = atom.AtomNameOnly()
z = 1 #atom.AtomicNumber() # Using 1 as default, might be updated for specific calculations
q = atom.Charge()
pos = atom.Position()
xc, yc, zc = cry.FractionalToCartesian(*pos)
if q >= 0.0:
r = 0.7
color = 'red'
M = 1.0
else:
r = 1.4
color = 'blue'
M = 1.0
sites.append([name, label, z, M, q, r, color, pos])
inf = EWALD(sites, lattice_parameters, Rlatt, gij, Rgij, volume, Rvolume, ew_alpha, prec)
print("")
print(f"Time for real space sum : {inf['time_real_space']:6}")
print(f"Time for real reciprocal sum: {inf['time_reciprocal_space']:6}")
print(f"Total time : {inf['time_total']:6}")
print("")
print("Madelung potential (electrostatic potential):")
MP_tot = 0.0
Ke = e * e / 4.0 / pi / e0 # J m
for isite in range(len(sites)):
name, label, z, M, q, r, color, pos = sites[isite]
UC1 = inf['UC1_list'][isite] # Potential in V
UC2 = inf['UC2_list'][isite] # Potential in V
UC3 = inf['UC3_list'][isite] # Potential in V
MP = inf['MP_list'][isite] # Potential in V
# Total Madelung energy in Joules is sum of q_i * V_i
# Here, q is in elementary charge units, so multiply by 'e' (Coulombs) to get J
# Then MP_tot is total energy in J.
MP_tot += q * MP * e # q (elementary charge units) * MP (Volts) * e (C/elementary charge) = Joules
# Output MP in Joules and eV.
# Note: Potential in Volts is numerically equal to Potential in eV/e (eV per elementary charge).
# 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).
print(f"{name:4}: {label:6}: q={q:8.2g}: "
+ f"MP = {MP * e:12.6g} J (= {UC1 * e:12.6g} + {UC2 * e:12.6g} - {UC3 * e:12.6g})")
print(f" MP = {MP:12.6g} eV (= {UC1:12.6g} + {UC2:12.6g} - {UC3:12.6g})")
MP_tot *= 0.5 # Total Madelung energy calculation usually requires 1/2 sum(qi * Vi)
print(f"Total Madelung energy in unit cell: {MP_tot:12.6g} J")
print(f"Total Madelung energy in unit cell: {MP_tot / e:12.6g} eV") # Convert total J to eV
# Charge is represented by q0 to define Madelung constant
# Lattice parameter a is represented by q0 to define Madelung constant
namej, labelj, zj, Mj, qj, rj, colorj, pos_j = sites[0]
print("")
print("Madelung constant")
print(f"NOTE: The a-axis length is taken as the representative atomic distance: a = {lattice_parameters[0]}")
print(f" The charge of the 0-th ion is taken as the representative ion charge: q = {qj}")
print(f"NOTE: This value is in the unit cell chemical formula")
print(" The following value must be devided by Z to get the Madeluing constant in the standard definition")
# Madelung constant alpha = - (MP_total / N) * (4*pi*e0*a / q^2)
# Total Madelung energy = Sum(0.5 * q_i * V_i)
# V_i = (1/(4*pi*e0)) * alpha * q_i / a
# E_M = Sum(0.5 * q_i * (1/(4*pi*e0)) * alpha_i * q_i / a)
# In the code, `MP_tot` is the total energy in Joules.
# We want Madelung constant `alpha_M = - E_M * a / (Ke * N * q^2)`.
# `E_M` is `MP_tot` (total energy in J).
# `a` is `lattice_parameters[0] * 1.0e-10` (meters).
# `Ke` is `e^2 / (4*pi*e0)` (J m).
# `qj` is the charge number, so `qj^2` is the square of the charge number.
# The definition `Madelung constant in unit cell: {:14.8g}".format(-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10))`
# 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`.
# `alpha = - E_M * (4*pi*e0 * a) / (q_ref^2)`.
# `E_M` is `MP_tot` (J).
# `a` is `lattice_parameters[0] * 1.0e-10` (m).
# `q_ref` is `qj * e` (C).
# So `alpha = -MP_tot * (4*pi*e0 * a) / (qj * e)^2`.
# `Ke = e^2 / (4*pi*e0)`. So `1 / (4*pi*e0) = Ke / e^2`.
# `alpha = -MP_tot * a / (Ke * qj^2)`.
# The existing calculation is `-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10`.
# This becomes `-MP_tot / (e^2 / (4*pi*e0)) * e / abs(qj) * a`.
# Which simplifies to `-MP_tot * (4*pi*e0) / e / abs(qj) * a`.
# This is not `alpha = -MP_tot * (4*pi*e0 * a) / (qj * e)^2`.
# There might be a specific definition of Madelung constant they are using with `Ke` and `e`.
# Let's assume the existing line of code for Madelung constant calculation is correct based on internal conventions.
print("Madelung constant in unit cell: {:14.8g}".format(-MP_tot / Ke * e / abs(qj) * lattice_parameters[0] * 1.0e-10))
app.terminate(pause = True)
if __name__ == '__main__':
main()