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

xrd_pymatgen.py をダウンロード

xrd_pymatgen.py
xrd_pymatgen.py
  1#Pymatgenチュートリアル⑥ XRDのシミュレーションをする https://qiita.com/ojiya/items/1b154c3698cff91c8a2b
  2
  3"""
  4X線回折（XRD）パターンシミュレーションスクリプト。
  5
  6概要:
  7    Pymatgenライブラリを使用して結晶のX線回折パターンをシミュレーションし、結果をプロット・保存します。
  8
  9詳細説明:
 10    本スクリプトは、CIFファイルから結晶構造を読み込み、指定されたX線源（例: CuKa1）と2θ範囲でXRDパターンを計算します。
 11    計算された個々の回折ピークに対してガウス関数を適用し、指定された半値半幅 (FWHM) でスムージングされたスペクトルを生成します。
 12    結果の回折ピーク情報と生成されたスペクトルデータは、標準出力、指定されたログファイル、およびExcelファイルに保存されます。
 13    また、計算されたスペクトルと個々のピークはmatplotlibによってプロットされ、視覚的に表示されます。
 14    コマンドライン引数により、入力CIFファイル、X線源、2θの最小値・最大値・ステップ、FWHMなどを柔軟に設定できます。
 15
 16関連リンク:
 17    - Pymatgen XRDCalculator ドキュメント:
 18      https://pymatgen.org/pymatgen.analysis.diffraction.xrd.html
 19    - Pymatgenチュートリアル（Qiita）:
 20      https://qiita.com/ojiya/items/1b154c3698cff91c8a2b
 21    - Pymatgen XRDソースコード (GitHub):
 22      https://github.com/materialsproject/pymatgen/blob/v2023.5.10/pymatgen/analysis/diffraction/xrd.py
 23"""
 24
 25import math
 26from math import exp
 27import numpy as np
 28import scipy as sp
 29import scipy.special
 30import matplotlib.pyplot as plt
 31
 32from pymatgen.core.periodic_table import Element
 33from pymatgen.io.cif import CifParser
 34from pymatgen.core.structure import Structure
 35from pymatgen.core.lattice import Lattice
 36from pymatgen.analysis.diffraction.xrd import XRDCalculator, WAVELENGTHS
 37from pymatgen.analysis.diffraction.tem import TEMCalculator
 38from pymatgen.analysis.diffraction.neutron import NDCalculator
 39
 40from tklib.tkapplication import tkApplication
 41from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
 42from tklib.tkvariousdata import tkVariousData
 43from tklib.tkgraphic.tkplotevent import tkPlotEvent
 44
 45
 46infile = 'SrTiO3.cif'
 47
 48# URL: https://github.com/materialsproject/pymatgen/blob/v2023.5.10/pymatgen/analysis/diffraction/xrd.py
 49# 'CuKa', 'CuKa2', 'CuKa1', 'CuKb1', 'MoKa', 'MoKa2', 'MoKa1', 'MoKb1',
 50# 'CrKa', 'CrKa2', 'CrKa1', 'CrKb1', 'FeKa', 'FeKa2', 'FeKa1', 'FeKb1', 'CoKa', 'CoKa2', 'CoKa1', 'CoKb1',
 51# 'AgKa', 'AgKa2', 'AgKa1', 'AgKb1'
 52"""
 53WAVELENGTHS = {
 54    "CuKa": 1.54184,
 55    "CuKa2": 1.54439,
 56    "CuKa1": 1.54056,
 57    "CuKb1": 1.39222,
 58    "MoKa": 0.71073,
 59    "MoKa2": 0.71359,
 60    "MoKa1": 0.70930,
 61    "MoKb1": 0.63229,
 62    "CrKa": 2.29100,
 63    "CrKa2": 2.29361,
 64    "CrKa1": 2.28970,
 65    "CrKb1": 2.08487,
 66    "FeKa": 1.93735,
 67    "FeKa2": 1.93998,
 68    "FeKa1": 1.93604,
 69    "FeKb1": 1.75661,
 70    "CoKa": 1.79026,
 71    "CoKa2": 1.79285,
 72    "CoKa1": 1.78896,
 73    "CoKb1": 1.63079,
 74    "AgKa": 0.560885,
 75    "AgKa2": 0.563813,
 76    "AgKa1": 0.559421,
 77    "AgKb1": 0.497082,
 78}
 79"""
 80
 81Xray_source     = 'CuKa1'
 82
 83Q2min  = 10.0
 84Q2max  = 80.0
 85Q2step = 0.02
 86
 87fwhm = 0.05
 88
 89figsize = (8, 4)
 90fontsize = 12
 91
 92
 93app    = tkApplication()
 94infile      = getarg( 1, infile)
 95Xray_source = getarg( 2, Xray_source)
 96Q2min       = getfloatarg( 3, Q2min)
 97Q2max       = getfloatarg( 4, Q2max)
 98Q2step      = getfloatarg( 5, Q2step)
 99fwhm        = getfloatarg( 6, fwhm)
100
101
102#==========================================
103# Main prgram
104#==========================================
105def Gaussian(x, x0, whalf):
106    """
107    ガウス関数を計算します。
108
109    概要:
110        ピークの中心と半値半幅を用いてガウス分布の値を計算します。
111        この関数はXRDパターンシミュレーションにおいて、各回折ピークの形状を表現するために使用されます。
112
113    Parameters:
114        :param x: float or numpy.ndarray - ガウス関数を評価する点または点の配列。
115        :param x0: float - ガウス関数のピーク中心（位置）。
116        :param whalf: float - ガウス関数の半値半幅 (FWHM)。
117
118    Returns:
119        :returns: float or numpy.ndarray - `x`におけるガウス関数の値。
120    """
121    a = whalf / 0.832554611
122    X = (x - x0) / a
123    return exp(-X * X)
124
125
126#==========================================
127# Main routine
128#==========================================
129def main():
130    """
131    XRDパターンシミュレーションのメイン処理を実行します。
132
133    概要:
134        CIFファイルから結晶構造を読み込み、X線回折パターンを計算して可視化・保存します。
135
136    詳細説明:
137        1. コマンドライン引数またはデフォルト設定から、入力CIFファイル、X線源、
138           2θ範囲、スペクトルの半値半幅(FWHM)などのパラメータを読み込みます。
139        2. ログファイルおよび出力Excelファイルのパスを設定し、標準出力をログファイルにもリダイレクトします。
140        3. Pymatgenの`CifParser`でCIFファイルを解析し、結晶構造オブジェクトを取得します。
141        4. `XRDCalculator`を使用して、指定されたX線源と2θ範囲で回折ピークデータを計算します。
142        5. 各回折ピークに対して`Gaussian`関数を適用し、スムージングされたXRDスペクトルを生成します。
143        6. 計算された回折ピークの詳細と生成されたスペクトルデータを標準出力とExcelファイルに保存します。
144        7. matplotlibを使用して、スムージングされたXRDスペクトルと個々の回折ピークをプロットし、表示します。
145           プロットにはインタラクティブなイベントハンドラが設定されており、カーソル位置のデータ情報を表示できます。
146        8. スクリプトの終了時に一時停止します。
147
148    Parameters:
149        None (コマンドライン引数から設定を読み込みます)
150
151    Returns:
152        None
153    """
154    logfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-out.txt"])
155    print(f"Open logfile [{logfile}]")
156    app.redirect(targets = ["stdout", logfile], mode = 'w')
157
158    outxlsxfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-pxrd.xlsx"])
159    nQ2 = int((Q2max - Q2min) / Q2step + 1.000001)
160    dQ2 = 6.0 * fwhm
161    nconv = int(dQ2 / Q2step + 1.00001)
162
163    print("")
164    print(f"input: {infile}")   
165    print(f"log file: {logfile}")
166    print(f"(output xlsx file: {outxlsxfile}")
167    print(f"X-ray source: {Xray_source}")
168    try:
169        wl = float(Xray_source)
170    except:
171        wl = WAVELENGTHS.get(Xray_source, None)
172    if wl:
173        print(f"  wavelength: {wl} angstrom")
174    print(f"2Theta range: {Q2min:9.5g} - {Q2max:9.5g} at {Q2step:9.5g} degree ({nQ2} points)")
175    print(f"2Theta range for a peak: {dQ2:9.5g} degree for {nconv} points")
176    print(f"  FWHM of Gaussian functin: {fwhm} degree")
177
178    x = np.arange(10,60,0.02)
179
180    parser = CifParser(infile)
181    mat = parser.get_structures(primitive=False, symmetrized=False)[0]
182
183    xrd = XRDCalculator(wavelength = Xray_source) #, debye_waller_factors = {Element('Sr'): 0.2, Element('O'): 0.6})
184    Q2min_diffraction = max([  0.0, Q2min - dQ2])
185    Q2max_diffraction = min([180.0, Q2max + dQ2])
186    diffractions = xrd.get_pattern(mat, two_theta_range = (Q2min_diffraction, Q2max_diffraction))
187#diffractions.as_dict()
188
189    ndiffractions = len(diffractions.x)
190    Q2   = diffractions.x
191    dhkl = diffractions.d_hkls
192    hkl  = [diffractions.hkls[i][0]['hkl'] for i in range(ndiffractions)]
193    mul  = [diffractions.hkls[i][0]['multiplicity'] for i in range(ndiffractions)]
194# y includes mul, LP and T
195    Int  = [diffractions.y[i] for i in range(ndiffractions)]
196
197    print("")
198    print(f"Diffractions: {ndiffractions} diffractions")
199    print(f"{'h':2} {'k':2} {'l':2}   {'m':2}    {'dhkl':8} {'2Theta':8} {'Intensity'}")
200    for i in range(len(Q2)):
201        h = hkl[i][0]
202        k = hkl[i][1]
203        if len(hkl[i]) == 4:
204            l = hkl[i][3]
205        else:
206            l = hkl[i][2]
207        print(f"{h:2} {k:2} {l:2}   {mul[i]:2}    {dhkl[i]:8.6f} {Q2[i]:8.6f} {Int[i]:8.4f}")
208#        print("  diffractions.hkls[i]=", diffractions.hkls[i])
209
210    print("")
211    print("Calculate specrum")
212    xQ2     = np.arange(Q2min, Q2max + Q2step, Q2step)
213    xrd_cal = np.zeros(nQ2)
214    for idf in range(ndiffractions):
215        idx0 = int((Q2[idf] - Q2min) / Q2step)
216
217        for j in range(idx0 - nconv, idx0 + nconv + 1):
218            if j < 0 or j >= nQ2:
219                continue
220
221            xrd_cal[j] += Int[idf] * Gaussian(xQ2[j], Q2[idf], fwhm)
222
223    print("")
224    print(f"Save to [{outxlsxfile}]")
225    tkVariousData().to_excel(outxlsxfile, ['2Theta (degree)', 'Intensity'], [xQ2, xrd_cal])
226
227    mode = ''
228    if mode == 'xrd':
229        xrd_pattern = xrd.plot_structures(
230            [mat],
231            two_theta_range=[Q2min, Q2max],
232            show = True,
233#            savefig="si_xrd.pdf"
234        )
235        exit()
236    elif mode == 'nd':
237        nd = NDCalculator()
238        nd_pattern = nd.plot_structures(
239            [mat],
240            two_theta_range=[Q2min, Q2max],
241            show = True,
242#            savefig="si_nd.pdf"
243        )
244        exit()
245    elif mode == 'tem':
246        tem = TEMCalculator()
247        tem_pattern = tem.get_plot_2d(
248            mat
249            )
250        tem_pattern.write_image("si_tem.pdf")
251        exit()
252    
253#=====================
254# plot
255#=====================
256
257    fig = plt.figure(figsize = figsize)
258    plot_event = tkPlotEvent(plt)
259
260    ax = fig.add_subplot(111)
261    ax.set_title(f"{infile} - source: {Xray_source}")
262    ax.tick_params(labelsize = fontsize)
263
264    ax.plot(xQ2, xrd_cal, color = 'black', linewidth = 0.5)
265    data = ax.plot(Q2, Int, linestyle = '', marker = 'o', markerfacecolor = 'red', markeredgecolor = 'red', markersize = 3.0)
266
267#plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": [Q2, dhkl, hkl, mul, Int]})
268    hkl_str = [f"{hkl[i][0]} {hkl[i][1]} {hkl[i][2]}" for i in range(ndiffractions)]
269    plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": data,
270                                "x": Q2, 'y': Int, "xlabel": "2Theta", "ylabel": "Intensity",
271                                 "xlist": [Q2, dhkl,  hkl_str, mul, Int],
272                                 "xlabels": ['Q2', 'dhkl', 'hkl', 'multiplicity', 'Intensity']
273                                 })
274
275    ax.set_xlim(Q2min, Q2max)
276
277    ax.set_xlabel(r'2$\theta$', fontsize = fontsize)
278    ax.set_ylabel('Intensity', fontsize = fontsize)
279
280    plot_event.register_event(fig, event = "button_press_event", 
281                    callback = lambda event: plot_event.onclick(event))
282
283    plt.tight_layout()
284    plt.pause(0.1)
285
286    app.terminate(pause = True)
287
288
289if __name__ == "__main__":
290    main()