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

XRD_diffraction_angles.py をダウンロード

XRD_diffraction_angles.py
XRD_diffraction_angles.py
  1"""
  2X線回折のブラッグ角を計算するスクリプト。
  3
  4詳細説明:
  5このスクリプトは、CIF (Crystallographic Information File) 形式の結晶構造データから、
  6X線回折におけるブラッグ角 (2θ)、面間隔 (d)、およびミラー指数 (hkl) を計算します。
  7指定されたX線源の波長と最大2θ角に基づき、可能なすべての回折ピークを探索し、
  8結果をソートして標準出力およびログファイルに出力します。
  9
 10関連リンク: :doc:`XRD_diffraction_angles_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
 17from mpl_toolkits.mplot3d import Axes3D
 18import matplotlib.pyplot as plt
 19
 20
 21from tklib.tkapplication import tkApplication
 22from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
 23from tklib.tkfile import tkFile
 24from tklib.tkcrystal.tkcif import tkCIF
 25from tklib.tkcrystal.tkcrystal import tkCrystal
 26from tklib.tkcrystal.tkatomtype import tkAtomType
 27from tkcrystalbase import *
 28
 29
 30WAVELENGTHS = {
 31    "CuKa": 1.54184,
 32    "CuKa2": 1.54439,
 33    "CuKa1": 1.54056,
 34    "CuKb1": 1.39222,
 35    "MoKa": 0.71073,
 36    "MoKa2": 0.71359,
 37    "MoKa1": 0.70930,
 38    "MoKb1": 0.63229,
 39    "CrKa": 2.29100,
 40    "CrKa2": 2.29361,
 41    "CrKa1": 2.28970,
 42    "CrKb1": 2.08487,
 43    "FeKa": 1.93735,
 44    "FeKa2": 1.93998,
 45    "FeKa1": 1.93604,
 46    "FeKb1": 1.75661,
 47    "CoKa": 1.79026,
 48    "CoKa2": 1.79285,
 49    "CoKa1": 1.78896,
 50    "CoKb1": 1.63079,
 51    "AgKa": 0.560885,
 52    "AgKa2": 0.563813,
 53    "AgKa1": 0.559421,
 54    "AgKb1": 0.497082,
 55}
 56
 57
 58infile = 'SrTiO3.cif'
 59
 60Xray_source = 'CuKa1'
 61
 62# G min to remove the origin of reciprocal space
 63Gmin = 1.0e-5
 64
 65# 2Theta max
 66Q2max = 150.0       # degree in 2Theta
 67
 68# Structure
 69# Lattice parameters (angstrom and degree)
 70#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
 71#sites = [
 72#     ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red',  np.array([0.0, 0.0, 0.0])]
 73#    ]
 74
 75
 76figsize = (8, 4)
 77fontsize = 12
 78
 79
 80#==========================================
 81# Main prgram
 82#==========================================
 83app    = tkApplication()
 84infile      = getarg( 1, infile)
 85Xray_source = getarg( 2, Xray_source)
 86Q2max       = getfloatarg(3, Q2max)
 87
 88
 89def main():
 90    """
 91    プログラムの主要な処理を実行し、X線回折のブラッグ角を計算して結果を出力します。
 92
 93    詳細説明:
 94    この関数は、以下のステップでX線回折角の計算と出力を行います。
 95    1. コマンドライン引数から入力ファイル名、X線源、最大2θ角を取得します。
 96    2. ログファイルをオープンし、標準出力もそこにリダイレクトします。
 97    3. 入力CIFファイルを読み込み、結晶構造情報を取得します。
 98    4. 単位格子の情報（格子定数、格子系、格子ベクトルなど）を表示します。
 99    5. 結晶サイト情報（原子の種類、位置など）を処理します。
100    6. 指定された最大2θ角に基づき、回折可能なhkl範囲を決定します。
101    7. 各hkl指数について、逆格子空間での距離Gを計算し、それから面間隔dと2θ角を計算します。
102    8. 計算された回折角、面間隔、hkl指数をリストに格納し、2θ角でソートします。
103    9. 重複する2θ角を持つピークをまとめ、多重度とともに整形して出力します。
104    10. プログラム終了時に一時停止します。
105
106    :returns: None
107    """
108    global lattice_parameters, sites
109
110    logfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-out.txt"])
111    print(f"Open logfile [{logfile}]")
112    app.redirect(targets = ["stdout", logfile], mode = 'w')
113
114    print("")
115    print(f"input: {infile}")   
116    print(f"log file: {logfile}")
117    print(f"X-ray source: {Xray_source}")
118    try:
119        wl = float(Xray_source)
120    except:
121        wl = WAVELENGTHS.get(Xray_source, None)
122    if wl:
123        print(f"  wavelength: {wl} angstrom")
124    print(f"2Theta max: {Q2max:9.5g} degree")
125
126    print("")
127    print("Read [{}]".format(infile))
128    cif = tkCIF()
129    cif.debug = False
130    cifdata = cif.ReadCIF(infile, find_valid_structure = True)
131    cif.Close()
132    if not cifdata:
133        terminate("Error: Could not get cifdat from infile [{}]".format(infile))
134
135#        cifdata.print()
136    cry = cifdata.GetCrystal()
137    cry.PrintInf()
138
139    lattice_parameters = cry.LatticeParameters()
140    ls = cry.lattice_system()
141    la = cry.lattice_axis()
142    aij = cry.LatticeVectors()
143    gij, Rgij = cry.Metrics()
144    volume = cry.Volume()
145
146    Raij = cry.ReciprocalLatticeVectors()
147    Rlatt = cry.ReciprocalLatticeParameters()
148    Rvolume = cry.ReciprocalVolume()
149
150    print("")
151    print("Lattice parameters:", lattice_parameters)
152    print(f"Lattice system: {ls}")
153    print(f"Unitcell axis system: {la}")
154
155    """
156    print("Lattice vectors:")
157    print("  ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
158    print("  ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
159    print("  az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
160    print("Metric tensor:")
161    print("  gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
162    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
163    print("       ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
164    print("Unit cell volume: {:12.4g} A^3".format(volume))
165
166    print("")
167    print("Reciprocal lattice parameters:", Rlatt)
168    print("Reciprocal lattice vectors:")
169    print("  Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
170    print("  Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
171    print("  Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
172    print("Reciprocal lattice metric tensor:")
173    print("  Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
174    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
175    print("        ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
176    print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
177    """
178
179    sites = []
180    AtomTypes = cry.AtomTypeList()
181    AtomSites = cry.ExpandedAtomSiteList()
182    allsites = []
183    for atom in AtomSites:
184        label  = atom.Label()
185        name   = atom.AtomNameOnly()
186        z      = 1 #atom.AtomicNumber()
187        q      = atom.Charge()
188        pos    = atom.Position()
189        xc, yc, zc = cry.FractionalToCartesian(*pos)
190        if q >= 0.0:
191            r = 0.7
192            color = 'red'
193            M = 1.0
194        else:
195            r = 1.4
196            color = 'blue'
197            M = 1.0
198
199        sites.append([name, label, z, M, q, r, color, pos])
200
201    dmin = wl / 2.0 / sin(0.5 * Q2max * torad)
202    hmax = int(lattice_parameters[0] / dmin)
203    kmax = int(lattice_parameters[1] / dmin)
204    lmax = int(lattice_parameters[2] / dmin)
205
206    print("")
207    print("hkl range:", hmax, kmax, lmax)
208
209# Calculate diffraction angles and store them in qlist list variable
210    org = np.array([0.0, 0.0, 0.0])
211    qlist = []
212    for l in range(-lmax, lmax+1):
213      for k in range(-kmax, kmax+1):
214        for h in range(-hmax, hmax+1):
215# Calculate distance in reciprocal space between (0, 0, 0) and (h, k, l)
216            G = distance(org, np.array([h, k, l]), Rgij)
217            if G < Gmin:
218                continue
219
220# Calculate lattice space from G
221            d = 1.0 / G
222
223            sinQ = wl / 2.0 / d
224            if sinQ >= 1.0:
225                continue
226            
227# Calculate diffraction angle 2Theta
228            Q2 = 2.0 * todeg * arcsin(sinQ)
229            if Q2 > Q2max:
230                continue
231
232            qlist.append([Q2, d, h, k, l])
233#            print("  2Q={:12.4g}  d={:12.6g}  ({:3d} {:3d} {:3d})".format(Q2, d, h, k, l))
234
235# Sort rlist by 2Theta (x[0] priority)
236    qlist.sort(key = lambda x: (x[0], x[2], x[3], x[4]))
237
238    peak_identical = {}
239    for qinf in qlist:
240        Q2, d, h, k, l = qinf
241        
242        h2, k2, l2 = cry.normalize_hkl(h, k, l)
243        if la == 'hexagonal' or la == 'trigonal':
244            key = f"{Q2:10.6f}:{h2}_{k2}_{-h2-k2}_{l2}"
245        else:
246            key = f"{Q2:10.6f}:{h2}_{k2}_{l2}"
247
248        if peak_identical.get(key, None):
249            peak_identical[key].append({ "Q2": Q2, "h": h, "k": k, "l": l, "d": d})
250        else:
251            peak_identical[key] = [{ "Q2": Q2, "h": h, "k": k, "l": l, "d": d}]
252
253    print("")
254    print("Diffraction angle, d, h, k, l:")
255    if la == 'hexagonal':
256#    print(f"{'h':2} {'k':2} {'l':2}    {'dhkl':10} {'2Theta':8}")
257        print(f"{'h':>2} {'k':>2} {'i':>2} {'l':>2}   {'m':>2}    {'dhkl':^10} {'2Theta':^8}")
258    else:
259        print(f"{'h':>2} {'k':>2} {'l':>2}   {'m':>2}    {'dhkl':^10} {'2Theta':^8}")
260    for key in sorted(peak_identical.keys()):
261        p_key = peak_identical[key]
262        p = peak_identical[key][len(p_key)-1]
263        h  = p['h']
264        k  = p['k']
265        l  = p['l']
266        m = cry.multiplicity_hkl(h, k, l)   #len(peak_identical[key])
267        Q2 = p['Q2']
268        d  = p['d']
269        if la == 'hexagonal':
270#        print(f"{h:2} {k:2} {l:2}    {d:10.6f} {Q2:8.6f}")
271            print(f"{h:2} {k:2} {-h-k:2} {l:2}   {m:2}    {d:10.6f} {Q2:8.6f}")
272        else:
273            print(f"{h:2} {k:2} {l:2}   {m:2}    {d:10.6f} {Q2:8.6f}")
274
275    print("")
276    app.terminate(pause = True)
277
278
279if __name__ == '__main__':
280    main()