diff --git a/.github/workflows/test-ubuntu-mamba.yml b/.github/workflows/test-ubuntu.yml similarity index 79% rename from .github/workflows/test-ubuntu-mamba.yml rename to .github/workflows/test-ubuntu.yml index 10e8d3eb2..e024c0130 100644 --- a/.github/workflows/test-ubuntu-mamba.yml +++ b/.github/workflows/test-ubuntu.yml @@ -1,4 +1,4 @@ -name: Test with Ubuntu, Mamba +name: Test with Ubuntu, Conda on: [push] @@ -19,7 +19,6 @@ jobs: uses: conda-incubator/setup-miniconda@v2.2.0 with: miniconda-version: latest - mamba-version: "*" channels: conda-forge,defaults channel-priority: true activate-environment: anaconda-client-env @@ -27,7 +26,7 @@ jobs: - name: Install conda dependencies run: | printenv - mamba install -y -c conda-forge python=${{ matrix.python-version }} "numpy>=1.23" "scipy>=1.8" "matplotlib>=3.6" "h5py>=3.5" "wxpython>=4.2" scikit-image scikit-learn pycifrw pandas jupyter plotly pyparsing pytest pytest-cov coverage + conda install -y -c conda-forge python=${{ matrix.python-version }} "numpy>=1.23" "scipy>=1.8" "matplotlib>=3.6" "h5py>=3.5" "wxpython>=4.2" scikit-image scikit-learn pandas jupyter plotly pyparsing pytest pytest-cov coverage - name: Install xraylarch and other dependencies with pip run: | pip install ".[all]" diff --git a/.github/workflows/test-windows-mamba.yml b/.github/workflows/test-windows-mamba.yml deleted file mode 100644 index 47d3a7f0b..000000000 --- a/.github/workflows/test-windows-mamba.yml +++ /dev/null @@ -1,35 +0,0 @@ -name: Test with Windows, Mamba - -on: [push] - -jobs: - build: - runs-on: windows-latest - strategy: - max-parallel: 5 - fail-fast: false - matrix: - python-version: ['3.9', '3.10', '3.11', '3.12'] - steps: - - uses: actions/checkout@v4 - - name: Set up Python from Miniconda - uses: conda-incubator/setup-miniconda@v2.2.0 - with: - miniconda-version: latest - mamba-version: "*" - channels: conda-forge,defaults - channel-priority: true - activate-environment: anaconda-client-env - python-version: ${{ matrix.python-version }} - - name: Install conda dependencies - run: | - printenv - mamba install -y -c conda-forge python=${{ matrix.python-version }} "numpy>=1.23" "scipy>=1.8" "matplotlib>=3.6" "h5py>=3.5" "wxpython>=4.2" scikit-image scikit-learn pycifrw pandas jupyter plotly pyparsing pytest pytest-cov coverage - - name: Install xraylarch and other dependencies with pip - run: | - printenv - pip install ".[all]" - - name: Run test suite - run: | - cd tests - python -m pytest diff --git a/.github/workflows/test-windows.yml b/.github/workflows/test-windows.yml index 4586f8b0f..9678fb779 100644 --- a/.github/workflows/test-windows.yml +++ b/.github/workflows/test-windows.yml @@ -23,7 +23,7 @@ jobs: - name: Install conda dependencies run: | printenv - conda install -y -c conda-forge python=${{ matrix.python-version }} "numpy>=1.23" "scipy>=1.8" "matplotlib>=3.6" "h5py>=3.5" "wxpython>=4.2" scikit-image scikit-learn pycifrw pandas jupyter plotly pyparsing pytest pytest-cov coverage + conda install -y -c conda-forge python=${{ matrix.python-version }} "numpy>=1.23" "scipy>=1.8" "matplotlib>=3.6" "h5py>=3.5" "wxpython>=4.2" scikit-image scikit-learn pandas jupyter plotly pyparsing pytest pytest-cov coverage - name: Install xraylarch and other dependencies with pip run: | printenv diff --git a/README.rst b/README.rst index 0f539f82c..f3e29f61e 100644 --- a/README.rst +++ b/README.rst @@ -1,10 +1,8 @@ Larch: Data Analysis Tools for X-ray Spectroscopy and More ============================================================ -.. image:: https://github.com/xraypy/xraylarch/actions/workflows/test-ubuntu-mamba.yml/badge.svg - :target: https://github.com/xraypy/xraylarch/actions/workflows/test-ubuntu-mamba.yml -.. image:: https://github.com/xraypy/xraylarch/actions/workflows/test-windows-mamba.yml/badge.svg - :target: https://github.com/xraypy/xraylarch/actions/workflows/test-ubuntu-mamba.yml +.. image:: https://github.com/xraypy/xraylarch/actions/workflows/test-ubuntu.yml/badge.svg + :target: https://github.com/xraypy/xraylarch/actions/workflows/test-ubuntu.yml .. image:: https://github.com/xraypy/xraylarch/actions/workflows/test-windows.yml/badge.svg :target: https://github.com/xraypy/xraylarch/actions/workflows/test-windows.yml diff --git a/installers/GetLarch.bat b/installers/GetLarch.bat index 50792df35..fb82a29e9 100644 --- a/installers/GetLarch.bat +++ b/installers/GetLarch.bat @@ -19,7 +19,7 @@ echo ## basic mamba installed, running updates set PATH=%prefix%;%prefix%\bin;%prefix%\condabin;%prefix%\Scripts;%PATH% echo ## Installing basic python scipy packages -call %prefix%\Scripts\mamba install -yc conda-forge python==3.11.5 numpy scipy matplotlib h5py scikit-image scikit-learn pycifrw pandas jupyter plotly wxpython fabio pyfai pymatgen mkl_fft tomopy +call %prefix%\Scripts\mamba install -yc conda-forge python==3.11.5 numpy scipy matplotlib h5py scikit-image scikit-learn pandas jupyter plotly wxpython fabio pyfai pymatgen mkl_fft tomopy echo ## Installing xraylarch and dependencies from PyPI call %prefix%\Scripts\pip install xraylarch[larix] diff --git a/installers/GetLarch.sh b/installers/GetLarch.sh index 126e65c50..7d8ef2efe 100644 --- a/installers/GetLarch.sh +++ b/installers/GetLarch.sh @@ -20,7 +20,7 @@ condafile="Mambaforge-$uname-x86_64.sh" logfile=GetLarch.log ## set list of conda packages to install from conda-forge -cforge_pkgs="python=>3.12.3 numpy scipy matplotlib h5py scikit-image scikit-learn pycifrw pandas jupyter plotly wxpython fabio pyfai pymatgen mkl_fft tomopy" +cforge_pkgs="python=>3.12.3 numpy scipy matplotlib h5py scikit-image scikit-learn pandas jupyter plotly wxpython fabio pyfai pymatgen mkl_fft tomopy" unset CONDA_EXE CONDA_PYTHON_EXE CONDA_PREFIX PROJ_LIB diff --git a/installers/conda_constructor/construct.yaml b/installers/conda_constructor/construct.yaml index 75ed52e67..f687b0542 100644 --- a/installers/conda_constructor/construct.yaml +++ b/installers/conda_constructor/construct.yaml @@ -54,5 +54,4 @@ specs: - wxpython>=4.2.0 - tomopy - pymatgen - - pycifrw - numdifftools diff --git a/larch/interpreter.py b/larch/interpreter.py index 8ca991bde..ff7dcd5b7 100644 --- a/larch/interpreter.py +++ b/larch/interpreter.py @@ -14,7 +14,8 @@ import numpy from pathlib import Path from copy import deepcopy - +from functools import partial +from inspect import signature from . import site_config from asteval import valid_symbol_name from .symboltable import SymbolTable, Group, isgroup @@ -154,9 +155,8 @@ def __init__(self, symtable=None, input=None, writer=None, group = self.symtable.set_symbol(groupname, value=Group(__name__=groupname)) for fname, fcn in list(entries.items()): - if callable(fcn): - setattr(group, fname, - Closure(func=fcn, _larch=self, _name=fname)) + if callable(fcn) and '_larch' in signature(fcn).parameters: + setattr(group, fname, Closure(func=fcn, _larch=self)) else: setattr(group, fname, fcn) diff --git a/larch/io/athena_project.py b/larch/io/athena_project.py index ec6b0f40b..3030bcdd1 100644 --- a/larch/io/athena_project.py +++ b/larch/io/athena_project.py @@ -18,6 +18,7 @@ from larch import Group, repr_value from larch import __version__ as larch_version from larch.utils import bytes2str, str2bytes, fix_varname, asfloat, unixpath +from larch.math import remove_dups from xraydb import guess_edge import asteval @@ -760,6 +761,8 @@ def read(self, filename=None, match=None, do_preedge=True, do_bkg=False, raise IOError("%s '%s': cannot find file" % (ERR_MSG, self.filename)) from larch.xafs import pre_edge, autobk, xftf + from larch.xafs.xafsutils import TINY_ENERGY + if not Path(filename).exists(): raise IOError("file '%s' not found" % filename) @@ -797,6 +800,17 @@ def read(self, filename=None, match=None, do_preedge=True, do_bkg=False, if not fnmatch(gname.lower(), match): continue this = getattr(data, gname) + this.energy = remove_dups(this.energy, tiny=TINY_ENERGY) + eorder = np.argsort(this.energy) + has_nan = False + for aname in dir(this): + obj = getattr(this, aname) + if isinstance(obj, np.ndarray) and len(eorder) == len(obj): + setattr(this, aname, obj[eorder]) + has_nan = has_nan or any(np.isnan(obj)) + this.energy = remove_dups(this.energy, tiny=TINY_ENERGY) + if has_nan: + print("NOTE: nans seen!!") this.athena_id = this.athena_params.id if use_hashkey: @@ -835,15 +849,18 @@ def read(self, filename=None, match=None, do_preedge=True, do_bkg=False, if is_chi: this.k = this.energy*1.0 this.chi = this.mu*1.0 + this.xdat = this.energy*1.0 + this.ydat = this.mu*1.0 del this.energy del this.mu + else: + this.xdat = 1.0*this.energy + this.ydat = 1.0*this.mu # add a selection flag and XAS datatypes, as used by Larix this.sel = 1 this.datatype = 'xas' this.filename = getattr(this, 'label', 'unknown') - this.xdat = 1.0*this.energy - this.ydat = 1.0*this.mu this.yerr = 1.0 this.plot_xlabel = 'energy' this.plot_ylabel = 'mu' diff --git a/larch/io/mergegroups.py b/larch/io/mergegroups.py index ac7021a40..2fb49a675 100644 --- a/larch/io/mergegroups.py +++ b/larch/io/mergegroups.py @@ -5,7 +5,7 @@ import os import numpy as np from larch import Group -from larch.math import interp, interp1d, index_of, remove_dups +from larch.math import interp, interp1d, index_of, remove_dups, remove_nans from larch.utils.logging import getLogger _logger = getLogger("larch.io.mergegroups") @@ -32,14 +32,20 @@ def merge_groups(grouplist, master=None, xarray='energy', yarray='mu', if master is None: master = grouplist[0] - xout = remove_dups(getattr(master, xarray)) + xarr = getattr(master, xarray) + dxperc = np.percentile(abs(np.diff(xarr)), [1, 2, 25]) + dxmin = min(0.01*dxperc[1], 1.e-4*dxperc[2]) + + + xout = remove_dups(xarr, dxmin) + xout = remove_nans(xout, interp=True) xmins = [min(xout)] xmaxs = [max(xout)] yvals = [] for g in grouplist: - x = getattr(g, xarray) - y = getattr(g, yarray) + x = remove_nans(remove_dups(getattr(g, xarray), dxmin), interp=True) + y = remove_nans(getattr(g, yarray), interp=True) yvals.append(interp(x, y, xout, kind=kind)) xmins.append(min(x)) xmaxs.append(max(x)) diff --git a/larch/wxlib/cif_browser.py b/larch/wxlib/cif_browser.py index ef057448c..a8a9c2eb7 100644 --- a/larch/wxlib/cif_browser.py +++ b/larch/wxlib/cif_browser.py @@ -23,9 +23,10 @@ from xraydb import atomic_number import larch + from larch import Group from larch.xafs import feff8l, feff6l -from larch.xrd.cif2feff import cif_sites +from larch.xrd.cif2feff import cif2feffinp, cif_cluster, site_label, fcompact from larch.utils import read_textfile, mkdir from larch.utils.paths import unixpath from larch.utils.strutils import fix_filename, unique_name, strict_ascii @@ -227,8 +228,8 @@ def createMainPanel(self): wids['feffvers'] = Choice(panel, choices=['6', '8'], default=1, size=(80, -1), action=self.onGetFeff) - wids['feff_site'] = Choice(panel, choices=['1', '2', '3', '4'], - size=(80, -1), + wids['feff_site'] = Choice(panel, choices=[' ', ' ', ' ', ' '], + size=(200, -1), action=self.onGetFeff) wids['feff_cluster_size'] = FloatSpin(panel, value=7.0, digits=2, increment=0.1, max_val=10, @@ -244,10 +245,10 @@ def createMainPanel(self): sizer.Add(SimpleText(panel, ' Absorbing Atom: '), (ir, 0), (1, 1), LEFT, 3) sizer.Add(wids['feff_central_atom'], (ir, 1), (1, 1), LEFT, 3) - sizer.Add(SimpleText(panel, ' Crystal Site: '), (ir, 2), (1, 1), LEFT, 3) - sizer.Add(wids['feff_site'], (ir, 3), (1, 1), LEFT, 3) - sizer.Add(SimpleText(panel, ' Edge: '), (ir, 4), (1, 1), LEFT, 3) - sizer.Add(wids['feff_edge'], (ir, 5), (1, 1), LEFT, 3) + sizer.Add(SimpleText(panel, ' Edge: '), (ir, 2), (1, 1), LEFT, 3) + sizer.Add(wids['feff_edge'], (ir, 3), (1, 1), LEFT, 3) + sizer.Add(SimpleText(panel, ' Crystal Site: '), (ir, 4), (1, 1), LEFT, 3) + sizer.Add(wids['feff_site'], (ir, 5), (1, 1), LEFT, 3) ir += 1 sizer.Add(SimpleText(panel, ' Cluster Size (\u212B): '), (ir, 0), (1, 1), LEFT, 3) @@ -292,16 +293,17 @@ def _swallow_plot_messages(s, panel=0): self.plotpanel.SetMaxSize((675, 400)) self.plotpanel.onPanelExposed = self.showXRD1D + fw_font = wx.Font(FONTSIZE, wx.MODERN, wx.NORMAL, wx.BOLD) + cif_panel = wx.Panel(rightpanel) wids['cif_text'] = wx.TextCtrl(cif_panel, value='', style=wx.TE_MULTILINE|wx.TE_READONLY, - size=(700, 450)) - wids['cif_text'].SetFont(Font(FONTSIZE+1)) + size=(725, 450)) + wids['cif_text'].SetFont(fw_font) cif_sizer = wx.BoxSizer(wx.VERTICAL) cif_sizer.Add(wids['cif_text'], 0, LEFT, 1) pack(cif_panel, cif_sizer) - self.nbpages = [] for label, page in (('CIF Text', cif_panel), ('1-D XRD Pattern', self.plotpanel), @@ -318,11 +320,12 @@ def _swallow_plot_messages(s, panel=0): wids['feff_text'] = wx.TextCtrl(feffinp_panel, value='', style=wx.TE_MULTILINE, - size=(700, 450)) + size=(725, 450)) + wids['feff_text'].SetFont(fw_font) wids['feff_text'].CanCopy() feffinp_panel.onPanelExposed = self.onGetFeff - wids['feff_text'].SetFont(Font(FONTSIZE+1)) + # wids['feff_text'].SetFont(Font(FONTSIZE+1)) feff_sizer = wx.BoxSizer(wx.VERTICAL) feff_sizer.Add(wids['feff_text'], 0, LEFT, 1) pack(feffinp_panel, feff_sizer) @@ -331,9 +334,10 @@ def _swallow_plot_messages(s, panel=0): wids['feffout_text'] = wx.TextCtrl(feffout_panel, value='', style=wx.TE_MULTILINE, - size=(700, 450)) + size=(725, 450)) wids['feffout_text'].CanCopy() - wids['feffout_text'].SetFont(Font(FONTSIZE+1)) + wids['feffout_text'].SetFont(fw_font) + # wids['feffout_text'].SetFont(Font(FONTSIZE+1)) feffout_sizer = wx.BoxSizer(wx.VERTICAL) feffout_sizer.Add(wids['feffout_text'], 0, LEFT, 1) pack(feffout_panel, feffout_sizer) @@ -448,13 +452,14 @@ def onShowCIF(self, event=None, cif_id=None): self.wids['feff_central_atom'].AppendItems(list(elems.keys())) self.wids['feff_central_atom'].Select(0) - el0 = list(elems.keys())[0] - edge_val = 'K' if atomic_number(el0) < 60 else 'L3' + catom = list(elems.keys())[0] + edge_val = 'K' if atomic_number(catom) < 60 else 'L3' self.wids['feff_edge'].SetStringSelection(edge_val) - sites = cif_sites(cif.ciftext, absorber=el0) + cluster = cif_cluster(ciftext=cif.ciftext, absorber=catom) + self.absorber_sites = cluster.all_sites[catom] try: - sites = ['%d' % (i+1) for i in range(len(sites))] + sites = list(self.absorber_sites.keys()) except: title = "Could not make sense of atomic sites" message = [f"Elements: {list(elems.keys())}", @@ -482,11 +487,14 @@ def onFeffCentralAtom(self, event=None): return catom = event.GetString() try: - sites = cif_sites(cif.ciftext, absorber=catom) - sites = ['%d' % (i+1) for i in range(len(sites))] + cluster = cif_cluster(ciftext=cif.ciftext, absorber=catom) + self.absorber_sites = cluster.all_sites[catom] + + sites = list(self.absorber_sites.keys()) self.wids['feff_site'].Clear() self.wids['feff_site'].AppendItems(sites) self.wids['feff_site'].Select(0) + except: self.write_message(f"could not get sites for central atom '{catom}'") title = f"Could not get sites for central atom '{catom}'" @@ -505,16 +513,37 @@ def onGetFeff(self, event=None): edge = self.wids['feff_edge'].GetStringSelection() version8 = '8' == self.wids['feffvers'].GetStringSelection() catom = self.wids['feff_central_atom'].GetStringSelection() - asite = int(self.wids['feff_site'].GetStringSelection()) + label = self.wids['feff_site'].GetStringSelection() + site_index = self.absorber_sites[label] csize = self.wids['feff_cluster_size'].GetValue() with_h = not self.wids['feff_without_h'].IsChecked() mineral = cif.get_mineralname() - folder = f'{catom:s}{asite:d}_{edge:s}_{mineral}_cif{cif.ams_id:d}' + + folder = f'{catom:s}{site_index:d}_{edge:s}_{mineral}_cif{cif.ams_id:d}' folder = unique_name(fix_filename(folder), self.feffruns_list) - fefftext = cif.get_feffinp(catom, edge=edge, cluster_size=csize, - absorber_site=asite, version8=version8, - with_h=with_h) + # extra titles from CIF + etitles = [f'CIF Source: Am Min Crystal Structure Database ID {cif.ams_id}', + f'Mineral Name: {cif.mineral.name.lower()}'] + pub = getattr(cif, 'publication', None) + if pub is not None: + authors = ', '.join(pub.authors) + ptext = f'{pub.journalname} {pub.volume} [{pub.page_first}:{pub.page_last}] ({pub.year})' + etitles.append(f'Publication: {ptext}') + etitles.append(f'Authors: {authors}') + + cell = [f'a={fcompact(cif.a)}', f'b={fcompact(cif.b)}', f'c={fcompact(cif.c)}', + f'alpha={fcompact(cif.alpha)}', f'beta={fcompact(cif.beta)}', + f'gamma={fcompact(cif.gamma)}'] + + etitles.append(f'Cell Parameters (A, degrees): {", ".join(cell)}') + etitles.append(f'Cell Volume (A^3): {cif.cell_volume:.5f}') + etitles.append(f'Crystal Density (gr/cm^3): {cif.crystal_density:.5f}') + etitles.append(f'Compound: {cif.compound}') + + fefftext = cif2feffinp(cif.ciftext, catom, edge=edge, cluster_size=csize, + absorber_site=site_index, version8=version8, + with_h=with_h, extra_titles=etitles) self.wids['feff_runfolder'].SetValue(folder) self.wids['feff_text'].SetValue(fefftext) @@ -570,7 +599,6 @@ def onRunFeff(self, event=None): wx.CallAfter(self.run_feff, folder, version8=version8) def run_feff(self, folder=None, version8=True): - # print("RUN FEFF ", folder) dname = Path(folder).name prog, cmd = feff8l, 'feff8l' if not version8: diff --git a/larch/wxxas/config.py b/larch/wxxas/config.py index 072c15efa..bf5c23b54 100644 --- a/larch/wxxas/config.py +++ b/larch/wxxas/config.py @@ -69,6 +69,7 @@ def make_array_choice(opts): NNORM_STRINGS = {int(v): k for k, v in NNORM_CHOICES.items()} NORM_METHODS = ('polynomial', 'mback') +PREEDGE_FORMS = {'Constant': 0, 'Linear': 1} ATHENA_CLAMPNAMES = {'none': 0, 'slight': 1, 'weak': 5, 'medium': 20, 'strong': 100, 'rigid': 500} @@ -88,7 +89,7 @@ def make_array_choice(opts): LARIX_PANELS = { 'xydata': - AnalysisTab('XY Data', 'larch.wxxas.xydata_panel.XYDataPanel', + AnalysisTab('XY Data', 'larch.wxxas.xydata_panel.XYDataPanel', 'Read and Manipulate XY Data from Column Data files'), 'xasnorm': AnalysisTab('XAS Normalization', 'larch.wxxas.xasnorm_panel.XASNormPanel', @@ -97,32 +98,32 @@ def make_array_choice(opts): AnalysisTab('Pre-edge Peaks', 'larch.wxxas.prepeak_panel.PrePeakPanel', 'Curve Fitting for XANES Pre-edge Peaks'), 'pca': - AnalysisTab('XAS PCA', 'larch.wxxas.pca_panel.PCAPanel', + AnalysisTab('XAS PCA', 'larch.wxxas.pca_panel.PCAPanel', 'Principal Component Analysis for XANES and EXAFS'), - 'lincombo': - AnalysisTab('XAS Linear Combo', 'larch.wxxas.lincombo_panel.LinearComboPanel', + 'lincombo': + AnalysisTab('XAS Linear Combo', 'larch.wxxas.lincombo_panel.LinearComboPanel', 'Linear Combination Analysis for XANES and EXAFS'), - 'regression': + 'regression': AnalysisTab('XAS Regression', 'larch.wxxas.regress_panel.RegressionPanel', 'Linear Regression and Feature Selection for XANES and EXAFS'), - 'exafs': - AnalysisTab('EXAFS', 'larch.wxxas.exafs_panel.EXAFSPanel', + 'exafs': + AnalysisTab('EXAFS', 'larch.wxxas.exafs_panel.EXAFSPanel', 'EXAFS Background Subtraction and Fourier Transforms'), - 'feffit': + 'feffit': AnalysisTab('FEFF Fitting', 'larch.wxxas.feffit_panel.FeffitPanel', 'EXAFS Path Fitting with FEFF calculations'), } LARIX_MODES = { 'all': ('All', [k for k in LARIX_PANELS]), - 'xydata': ('General XY Data Visualization and Fitting', ('xydata', 'lmfit')), + 'xydata': ('General XY Data Visualization and Fitting', ('xydata', 'lmfit')), 'xas': ('XANES and EXAFS', ('xasnorm', 'prepeaks', 'pca', 'lincombo', 'exafs', 'feffit')), 'exafs': ('EXAFS only', ('xasnorm', 'exafs', 'feffit')), 'xanes': ('XANES only', ('xasnorm', 'prepeaks', 'pca', 'lincombo')), 'xrf': ('XRF Mapping and Analysis', ('maproi', 'mapareas', 'maptomo', 'mapxrf')), 'xrd1d': ('XRD 1D', ('xrd1d', )), } - + class CVar: """configuration variable""" def __init__(self, name, value, dtype, choices=None, desc='a variable', @@ -273,6 +274,7 @@ def __repr__(self): CVar('pre1', -200, 'float', step=5, desc='low-energy fit range for pre-edge line,\nrelative to E0'), CVar('pre2', -30, 'float', step=5, desc='high-energy fit range for pre-edge line,\nrelative to E0'), CVar('nvict', 0, 'int', min=0, max=3, desc='Victoreen order for pre-edge fitting\n(Energy^(-nvict))'), + CVar('npre', 1, 'int', min=0, max=1, desc='Pre-edge form, 0 for Constant, 1 for Linear'), CVar('show_pre', False, 'bool', desc='whether to show pre-edge energy range (pre1, pre2)'), CVar('norm_method', 'polynomial', 'choice', choices=NORM_METHODS, desc='normalization method'), CVar('nnorm', 'linear', 'choice', choices=list(NNORM_CHOICES.keys()), diff --git a/larch/wxxas/feffit_panel.py b/larch/wxxas/feffit_panel.py index 4a4f389ae..ecc78de99 100644 --- a/larch/wxxas/feffit_panel.py +++ b/larch/wxxas/feffit_panel.py @@ -161,11 +161,12 @@ {groupname:s}.feffit_history.insert(0, _feffit_result) """ -COMMANDS['path2chi'] = """# generate chi(k) and chi(R) for each path +COMMANDS['path2chi'] = """# generate chi(k), chi(R), and chi(q) for each path for label, path in {paths_name:s}.items(): path.calc_chi_from_params({pargroup_name:s}) xftf(path, kmin={kmin:.3f}, kmax={kmax:.3f}, dk={dk:.3f}, window='{kwindow:s}', kweight={kweight:.3f}) + xftr(path, rmin={rmin:.3f}, rmax={rmax:.3f}, dr={dr:.3f}, window='{rwindow:s}') #endfor """ @@ -1041,13 +1042,24 @@ def process(self, dgroup=None, **kws): has_data = 1.e-12 < (tchi**2).sum() except: has_data = False - cmds = [COMMANDS['feffit_trans'].format(**conf)] + + cmds = [] _feffit_dataset = getattr(self.larch.symtable, '_feffit_dataset', None) if _feffit_dataset is None: cmds.append(COMMANDS['feffit_dataset_init']) if has_data: + cmds.append("## SET DATA GROUP ") + ftargs = dict(kmin=opts['fit_kmin'], kmax=opts['fit_kmax'], dk=opts['fit_dk'], + kwindow=opts['fit_kwindow'], kweight=opts['plot_kw'], + rmin=opts['fit_rmin'], rmax=opts['fit_rmax'], + dr=opts.get('fit_dr', 0.1), rwindow='hanning') + + cmds.append(COMMANDS['xft'].format(groupname=dgroup.groupname, **ftargs)) cmds.append(f"_feffit_dataset.set_datagroup({dgroup.groupname})") cmds.append(f"_feffit_dataset.refine_bkg = {opts['refine_bkg']}") + cmds.append(COMMANDS['feffit_trans'].format(**conf)) + + self.larch.eval('\n'.join(cmds)) return opts @@ -1150,12 +1162,10 @@ def onPlot(self, evt=None, dataset_name='_feffit_dataset', dgroup = dataset.data else: dgroup = self.controller.get_group() - # print("no data? dgroup ", dataset.data, dgroup) if dgroup is not None: self.larch.eval(f"{dataset_name}.set_datagroup({dgroup.groupname})") dataset = getattr(self.larch.symtable, dataset_name, None) dgroup = dataset.data - # print("now data now dgroup ", dataset, getattr(dataset, 'data', None), dgroup) opts = self.process(dgroup) opts.update(**kws) @@ -1702,7 +1712,9 @@ def onFitModel(self, event=None, dgroup=None): dgroup.feffit_history[0].timestamp = time.strftime("%Y-%b-%d %H:%M") dgroup.feffit_history[0].label = label - fitlabels = [fhist.label for fhist in dgroup.feffit_history[1:]] + fitlabels = [] + for ix, fhist in enumerate(dgroup.feffit_history[1:]): + fitlabels.append(getattr(fhist, 'label', f'fit {1+ix}')) if label in fitlabels: count = 1 while label in fitlabels: @@ -2093,12 +2105,11 @@ def onPlot(self, event=None): return dset = result.datasets[0] dgroup = dset.data - if not hasattr(dset.data, 'rwin'): - dset._residual(result.params) - dset.save_outputs() + trans = dset.transform dset.prepare_fit(result.params) dset._residual(result.params) + dset.save_outputs() opts = self.feffit_panel.read_form(dgroup=dgroup) opts = {'build_fitmodel': False} diff --git a/larch/wxxas/xasgui.py b/larch/wxxas/xasgui.py index 6dbb3d1b6..58c8591f9 100644 --- a/larch/wxxas/xasgui.py +++ b/larch/wxxas/xasgui.py @@ -802,7 +802,7 @@ def createMenus(self): "Browse CIF Structure, run Feff", self.onCIFBrowse) MenuItem(self, feff_menu, "Generate Feff input from general structures, Run Feff", "Generate Feff input from general structures, run Feff", self.onStructureBrowse) - MenuItem(self, feff_menu, "Browse Feff Calculations", + MenuItem(self, feff_menu, "Browse Feff Calculations, Import Paths", "Browse Feff Calculations, Get Feff Paths", self.onFeffBrowse) self.menubar.Append(feff_menu, "Feff") diff --git a/larch/wxxas/xasnorm_panel.py b/larch/wxxas/xasnorm_panel.py index 3b95a41fa..7969fbab4 100644 --- a/larch/wxxas/xasnorm_panel.py +++ b/larch/wxxas/xasnorm_panel.py @@ -26,7 +26,7 @@ from larch.wxlib.plotter import last_cursor_pos from .xas_dialogs import EnergyUnitsDialog from .taskpanel import TaskPanel, autoset_fs_increment, update_confval -from .config import (make_array_choice, EDGES, ATSYMS, +from .config import (make_array_choice, EDGES, ATSYMS, PREEDGE_FORMS, NNORM_CHOICES, NNORM_STRINGS, NORM_METHODS) np.seterr(all='ignore') @@ -78,12 +78,12 @@ def build_display(self): plot_one = Button(trow, 'Plot Current Group', size=(175, -1), action=self.onPlotOne) - self.plotsel_op = Choice(trow, choices=list(PlotSel_Choices.keys()), + self.plotsel_op = Choice(trow, choices=list(PlotSel_Choices), action=self.onPlotSel, size=(300, -1)) - self.plotone_op = Choice(trow, choices=list(PlotOne_Choices.keys()), + self.plotone_op = Choice(trow, choices=list(PlotOne_Choices), action=self.onPlotOne, size=(300, -1)) - self.plot_erange = Choice(trow, choices=list(Plot_EnergyRanges.keys()), + self.plot_erange = Choice(trow, choices=list(Plot_EnergyRanges), action=self.onPlotEither, size=(175, -1)) opts = {'digits': 2, 'increment': 0.05, 'value': 0, 'size': (FSIZE, -1)} @@ -153,7 +153,7 @@ def build_display(self): # step rows nnorm_panel = wx.Panel(panel) - self.wids['nnorm'] = Choice(nnorm_panel, choices=list(NNORM_CHOICES.keys()), + self.wids['nnorm'] = Choice(nnorm_panel, choices=list(NNORM_CHOICES), size=(150, -1), action=self.onNNormChoice, default=2) self.wids['auto_nnorm'] = Check(nnorm_panel, default=True, label='auto?', @@ -168,8 +168,11 @@ def build_display(self): action=self.onEnergyRef, size=(300, -1)) self.wids['nvict'] = Choice(panel, choices=('0', '1', '2', '3'), - size=(100, -1), action=self.onNormMethod, + size=(150, -1), action=self.onNormMethod, default=0) + self.wids['npre'] = Choice(panel, choices=list(PREEDGE_FORMS), + size=(150, -1), action=self.onNormMethod) + self.wids['npre'].SetSelection(1) opts = {'size': (FSIZE, -1), 'digits': 2, 'increment': 5.0, 'action': self.onSet_Ranges, @@ -222,6 +225,9 @@ def build_display(self): use_auto = Button(panel, 'Use Default Settings', size=(200, -1), action=self.onResetNorm) + compare_norms = Button(panel, 'Compare Normalization Methods', size=(200, -1), + action=self.onCompareNorm) + def CopyBtn(name): return Button(panel, 'Copy', size=(60, -1), action=partial(self.onCopyParam, name)) @@ -267,6 +273,8 @@ def CopyBtn(name): panel.Add(pre_panel, dcol=2) panel.Add(CopyBtn('xas_pre'), dcol=1, style=RIGHT) + panel.Add(SimpleText(panel, 'Pre-edge Type:'), newrow=True) + panel.Add(self.wids['npre'], dcol=2) panel.Add(SimpleText(panel, 'Victoreen order:'), newrow=True) panel.Add(self.wids['nvict'], dcol=2) @@ -274,13 +282,14 @@ def CopyBtn(name): panel.Add(HLine(panel, size=(HLINEWID, 3)), dcol=4, newrow=True) add_text('Normalization : ') - panel.Add(self.wids['norm_method'], dcol=2) + panel.Add(self.wids['norm_method'], dcol=1) + panel.Add(compare_norms) panel.Add(CopyBtn('xas_norm'), dcol=1, style=RIGHT) add_text('Norm Energy range: ') panel.Add(nor_panel, dcol=2) panel.Add(SimpleText(panel, 'Polynomial Type:'), newrow=True) - panel.Add(nnorm_panel, dcol=3) + panel.Add(nnorm_panel, dcol=2) panel.Add(HLine(panel, size=(HLINEWID, 3)), dcol=4, newrow=True) panel.Add(self.wids['is_frozen'], newrow=True) @@ -300,8 +309,6 @@ def get_config(self, dgroup=None): dgroup = self.controller.get_group() if dgroup is None: return self.get_defaultconfig() - self.read_form() - defconf = self.get_defaultconfig() conf = getattr(dgroup.config, self.configname, defconf) @@ -348,15 +355,12 @@ def get_config(self, dgroup=None): if conf['energy_ref'] in (None, 'None'): conf['energy_ref'] = fname - conf['energy_shift'] = getattr(dgroup,'energy_shift', conf['energy_shift']) - if hasattr(dgroup, 'e0') and conf['atsym'] == '?': atsym, edge = guess_edge(dgroup.e0) conf['atsym'] = atsym conf['edge'] = edge - if hasattr(dgroup, 'mback_params'): conf['atsym'] = getattr(dgroup.mback_params, 'atsym', conf['atsym']) conf['edge'] = getattr(dgroup.mback_params, 'edge', conf['edge']) @@ -369,11 +373,11 @@ def fill_form(self, dgroup): opts = self.get_config(dgroup) self.skip_process = True if self.is_xasgroup(dgroup): - if self.plotone_op.GetCount() != len(PlotOne_Choices.keys()): - self.plotone_op.SetChoices(list(PlotOne_Choices.keys())) + if self.plotone_op.GetCount() != len(PlotOne_Choices): + self.plotone_op.SetChoices(list(PlotOne_Choices)) self.plotone_op.SetSelection(1) - if self.plotsel_op.GetCount() != len(PlotSel_Choices.keys()): - self.plotsel_op.SetChoices(list(PlotSel_Choices.keys())) + if self.plotsel_op.GetCount() != len(PlotSel_Choices): + self.plotsel_op.SetChoices(list(PlotSel_Choices)) self.plotsel_op.SetSelection(1) groupnames = list(self.controller.file_groups.keys()) @@ -402,6 +406,7 @@ def fill_form(self, dgroup): self.wids['energy_shift'].SetValue(opts['energy_shift']) self.wids['nvict'].SetStringSelection("%d" % opts['nvict']) + self.wids['npre'].SetSelection(opts['npre']) self.wids['show_e0'].SetValue(opts['show_e0']) self.wids['auto_e0'].SetValue(opts['auto_e0']) self.wids['auto_nnorm'].SetValue(opts.get('auto_nnorm', 0)) @@ -416,13 +421,12 @@ def fill_form(self, dgroup): self.wids['show_pre'].SetValue(opts['show_pre']) self.wids['show_norm'].SetValue(opts['show_norm']) - frozen = opts.get('is_frozen', False) frozen = getattr(dgroup, 'is_frozen', frozen) self.wids['is_frozen'].SetValue(frozen) self._set_frozen(frozen) - wx.CallAfter(self.unset_skip_process) + self.skip_process = False def set_nnorm_widget(self, nnorm=None): conf = self.get_config() @@ -440,7 +444,7 @@ def set_nnorm_widget(self, nnorm=None): def unset_skip_process(self): self.skip_process = False - def read_form(self): + def read_form(self, who='x'): "read form, return dict of values" form_opts = {} form_opts['e0'] = self.wids['e0'].GetValue() @@ -449,11 +453,11 @@ def read_form(self): val = self.wids[attr].GetValue() if val == 0: val = None form_opts[attr] = val - form_opts['energy_shift'] = self.wids['energy_shift'].GetValue() form_opts['nnorm'] = NNORM_CHOICES.get(self.wids['nnorm'].GetStringSelection(), 1) form_opts['nvict'] = int(self.wids['nvict'].GetStringSelection()) + form_opts['npre'] = PREEDGE_FORMS.get(self.wids['npre'].GetStringSelection(), 1) form_opts['plotone_op'] = self.plotone_op.GetStringSelection() form_opts['plotsel_op'] = self.plotsel_op.GetStringSelection() form_opts['plot_voff'] = self.wids['plot_voff'].GetValue() @@ -480,8 +484,10 @@ def onNormMethod(self, evt=None): nnorm = get_auto_nnorm(self.wids['norm1'].GetValue(), self.wids['norm2'].GetValue()) + npre_sel = self.wids['nvict'].GetStringSelection() + npre = 0 if npre_sel.lower().startswith('con') else 1 nvict = int(self.wids['nvict'].GetStringSelection()) - self.update_config({'norm_method': method, 'nnorm': nnorm, 'nvict': nvict}) + self.update_config({'norm_method': method, 'nnorm': nnorm, 'nvict': nvict, 'npre': npre}) if method.startswith('mback'): dgroup = self.controller.get_group() cur_elem = self.wids['atsym'].GetStringSelection() @@ -491,7 +497,7 @@ def onNormMethod(self, evt=None): self.wids['atsym'].SetStringSelection(atsym) self.update_config({'edge': edge, 'atsym': atsym}) time.sleep(0.002) - wx.CallAfter(self.onReprocess) + self.onReprocess() def _set_frozen(self, frozen): try: @@ -516,20 +522,19 @@ def onEnergyRef(self, evt=None): dgroup.energy_ref = eref self.update_config({'energy_ref': eref}, dgroup=dgroup) - def onPlotEither(self, evt=None): + def onPlotEither(self, evt=None, process=True): if self.last_plot_type == 'multi': self.onPlotSel(evt=evt) else: - self.onPlotOne(evt=evt) + self.onPlotOne(evt=evt, process=process) - def onPlotOne(self, evt=None): + def onPlotOne(self, evt=None, process=True): self.last_plot_type = 'one' - self.plot(self.controller.get_group()) - wx.CallAfter(self.controller.set_focus) + self.plot(self.controller.get_group(), process=process) + self.controller.set_focus() def onVoffset(self, evt=None): - time.sleep(0.002) - wx.CallAfter(self.onPlotSel) + self.onPlotSel() def onPlotSel(self, evt=None): newplot = True @@ -607,7 +612,7 @@ def onPlotSel(self, evt=None): set_zoomlimits(ppanel, zoom_limits) or ppanel.unzoom_all() ppanel.canvas.draw() - wx.CallAfter(self.controller.set_focus) + self.controller.set_focus() def onAuto_NNORM(self, evt=None): if evt.IsChecked(): @@ -616,7 +621,7 @@ def onAuto_NNORM(self, evt=None): self.set_nnorm_widget(nnorm) self.wids['auto_nnorm'].SetValue(0) time.sleep(0.001) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onResetNorm(self, evt=None): auto_nnorm = self.wids['auto_nnorm'].GetValue() @@ -631,10 +636,18 @@ def onResetNorm(self, evt=None): self.wids['auto_e0'].SetValue(1) self.wids['auto_e0'].SetValue(1) self.wids['nvict'].SetSelection(0) + self.wids['npre'].SetSelection(1) for attr in ('pre1', 'pre2', 'norm1', 'norm2'): self.wids[attr].SetValue(defaults[attr]) self.onReprocess() + def onCompareNorm(self, evt=None): + dgroup = self.controller.get_group() + self.ensure_xas_processed(dgroup, force_mback=True) + plot_yarrays = self.get_plot_arrays(dgroup, plot_choice='mback_poly') + self.plot(dgroup, plot_yarrays=plot_yarrays) + + def onCopyAuto(self, evt=None): opts = dict(pre1=0, pre2=0, nvict=0, norm1=0, norm2=0, norm_method='polynomial', nnorm=2, auto_e0=1, @@ -651,12 +664,10 @@ def onCopyAuto(self, evt=None): def onSaveConfigBtn(self, evt=None): conf = self.get_config() conf.update(self.read_form()) - # self.set_defaultconfig(conf) def onCopyParam(self, name=None, evt=None): conf = self.get_config() - form = self.read_form() - conf.update(form) + conf.update(self.read_form()) dgroup = self.controller.get_group() self.update_config(conf) self.fill_form(dgroup) @@ -700,24 +711,21 @@ def onAuto_XASE0(self, evt=None): dgroup = self.controller.get_group() find_e0(dgroup) self.update_config({'e0': dgroup.e0}) - time.sleep(0.002) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onSet_XASE0(self, evt=None, value=None): "handle setting auto e0 / show e0" auto_e0 = self.wids['auto_e0'].GetValue() self.update_config({'e0': self.wids['e0'].GetValue(), 'auto_e0':self.wids['auto_e0'].GetValue()}) - time.sleep(0.002) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onSet_XASE0Val(self, evt=None, value=None): "handle setting e0" self.wids['auto_e0'].SetValue(0) self.update_config({'e0': self.wids['e0'].GetValue(), 'auto_e0':self.wids['auto_e0'].GetValue()}) - time.sleep(0.002) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onSet_EnergyShift(self, evt=None, value=None): conf = self.get_config() @@ -733,7 +741,7 @@ def onSet_EnergyShift(self, evt=None, value=None): this.energy_shift = this.config.xasnorm['energy_shift'] = eshift self.stale_groups.append(this) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onSet_XASStep(self, evt=None, value=None): "handle setting edge step" @@ -744,7 +752,7 @@ def onSet_XASStep(self, evt=None, value=None): self.update_config({'edge_step': abs(edge_step), 'auto_step': False}) autoset_fs_increment(self.wids['step'], abs(edge_step)) time.sleep(0.01) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onSet_Ranges(self, evt=None, **kws): @@ -752,8 +760,7 @@ def onSet_Ranges(self, evt=None, **kws): for attr in ('pre1', 'pre2', 'norm1', 'norm2'): conf[attr] = self.wids[attr].GetValue() self.update_config(conf) - time.sleep(0.01) - wx.CallAfter(self.onReprocess) + self.onReprocess() def pin_callback(self, opt='__', xsel=None, relative_e0=True, **kws): """ @@ -772,8 +779,7 @@ def pin_callback(self, opt='__', xsel=None, relative_e0=True, **kws): self.wids['auto_e0'].SetValue(0) elif opt in ('pre1', 'pre2', 'norm1', 'norm2'): self.wids[opt].SetValue(xsel-e0) - time.sleep(0.01) - wx.CallAfter(self.onReprocess) + self.onReprocess() def onReprocess(self, evt=None, value=None, **kws): "handle request reprocess" @@ -785,14 +791,12 @@ def onReprocess(self, evt=None, value=None, **kws): return if not hasattr(dgroup.config, self.configname): return - form = self.read_form() self.process(dgroup=dgroup) if self.stale_groups is not None: for g in self.stale_groups: self.process(dgroup=g, force=True) self.stale_groups = None - self.onPlotEither() - + self.onPlotEither(process=False) def process(self, dgroup=None, force_mback=False, force=False, use_form=True, **kws): """ handle process (pre-edge/normalize) of XAS data from XAS form @@ -806,7 +810,7 @@ def process(self, dgroup=None, force_mback=False, force=False, use_form=True, ** self.skip_process = True conf = self.get_config(dgroup) - form = self.read_form() + form = self.read_form('process') if not use_form: form.update(self.get_defaultconfig()) @@ -882,19 +886,19 @@ def process(self, dgroup=None, force_mback=False, force=False, use_form=True, ** if not form['auto_step']: copts.append("step=%s" % gformat(float(edge_step))) - for attr in ('pre1', 'pre2', 'nvict', 'nnorm', 'norm1', 'norm2'): + for attr in ('pre1', 'pre2', 'nvict', 'npre', 'nnorm', 'norm1', 'norm2'): val = form[attr] if val is None or val == 'auto': val = 'None' - elif attr in ('nvict', 'nnorm'): + elif attr in ('nvict', 'nnorm', 'npre'): if val in NNORM_CHOICES: val = NNORM_CHOICES[val] val = int(val) else: val = f"{float(val):.2f}" copts.append(f"{attr}={val}") - # print("process PreEdge ", copts) - self.larch_eval("pre_edge(%s)" % (', '.join(copts))) + + self.larch_eval("pre_edge(%s) #a" % (', '.join(copts))) self.larch_eval("{group:s}.norm_poly = 1.0*{group:s}.norm".format(**form)) if not hasattr(dgroup, 'e0'): self.skip_process = False @@ -932,6 +936,8 @@ def process(self, dgroup=None, force_mback=False, force=False, use_form=True, ** else: norm_expr = """{group:s}.norm = 1.0*{group:s}.norm_{normmeth:s} {group:s}.norm *= {group:s}.edge_step_{normmeth:s}/{edge_step:.8f}""" + + self.larch_eval(norm_expr.format(**form)) if norm_method.startswith('area'): @@ -966,11 +972,11 @@ def process(self, dgroup=None, force_mback=False, force=False, use_form=True, ** conf['atsym'] = getattr(dgroup.mback_params, 'atsym') conf['edge'] = getattr(dgroup.mback_params, 'edge') self.update_config(conf, dgroup=dgroup) - # print("process updated conf ", dgroup, conf) - wx.CallAfter(self.unset_skip_process) + + self.skip_process = False - def get_plot_arrays(self, dgroup): + def get_plot_arrays(self, dgroup, plot_choice=None): lab = plotlabels.norm if dgroup is None: return @@ -981,7 +987,10 @@ def get_plot_arrays(self, dgroup): req_attrs = ['e0', 'norm', 'dmude', 'd2mude', 'pre_edge'] - pchoice = PlotOne_Choices.get(self.plotone_op.GetStringSelection(), 'norm') + if plot_choice in PlotOne_Choices.values(): + pchoice = plot_choice + else: + pchoice = PlotOne_Choices.get(self.plotone_op.GetStringSelection(), 'norm') if pchoice in ('mu', 'norm', 'i0', 'flat', 'dmude', 'd2mude'): lab = getattr(plotlabels, pchoice) @@ -1064,8 +1073,10 @@ def get_plot_arrays(self, dgroup): ival = min(len(y4e0)-1, index_of(dgroup.energy, dgroup.e0 + val)) dgroup.plot_extras.append(('marker', dgroup.e0+val, y4e0[ival], popts)) + return dgroup.plot_yarrays + def plot(self, dgroup, title=None, plot_yarrays=None, yoff=0, - delay_draw=True, multi=False, new=True, with_extras=True, **kws): + delay_draw=True, multi=False, new=True, with_extras=True, process=True, **kws): if self.skip_plotting: return @@ -1082,11 +1093,11 @@ def plot(self, dgroup, title=None, plot_yarrays=None, yoff=0, if groupname is None: return - self.ensure_xas_processed(dgroup, force_mback=True) - self.get_plot_arrays(dgroup) + if process: + self.ensure_xas_processed(dgroup, force_mback=True) - if plot_yarrays is None and hasattr(dgroup, 'plot_yarrays'): - plot_yarrays = dgroup.plot_yarrays + if plot_yarrays is None: + plot_yarrays = self.get_plot_arrays(dgroup) popts = self.controller.get_plot_conf() popts.update(kws) diff --git a/larch/xafs/feffit.py b/larch/xafs/feffit.py index c59b77845..3d2b6652a 100644 --- a/larch/xafs/feffit.py +++ b/larch/xafs/feffit.py @@ -113,7 +113,7 @@ def make_karrays(self, k=None, chi=None): self.k_ = self.kstep * arange(self.nfft, dtype='float64') self.r_ = self.rstep * arange(self.nfft, dtype='float64') - def _xafsft(self, chi, group=None, rmax_out=10, **kws): + def _xafsft(self, chi, group=None, rmax_out=10, with_chiq=False, **kws): "returns " for key, val in kws: if key == 'kw': @@ -121,26 +121,37 @@ def _xafsft(self, chi, group=None, rmax_out=10, **kws): setattr(self, key, val) self.make_karrays() - out = self.fftf(chi) + outr = self.fftf(chi) irmax = int(min(self.nfft/2, 1.01 + rmax_out/self.rstep)) group = set_xafsGroup(group, _larch=self._larch) r = self.rstep * arange(irmax) - mag = sqrt(out.real**2 + out.imag**2) + mag = sqrt(outr.real**2 + outr.imag**2) group.kwin = self.kwin[:len(chi)] group.r = r[:irmax] - group.chir = out[:irmax] + group.chir = outr[:irmax] group.chir_mag = mag[:irmax] - group.chir_pha = complex_phase(out[:irmax]) - group.chir_re = out.real[:irmax] - group.chir_im = out.imag[:irmax] + group.chir_pha = complex_phase(outr[:irmax]) + group.chir_re = outr.real[:irmax] + group.chir_im = outr.imag[:irmax] if self.rwin is None: xmin = max(self.rbkg, self.rmin) self.rwin = ftwindow(self.r_, xmin=xmin, xmax=self.rmax, dx=self.dr, dx2=self.dr2, window=self.rwindow) group.rwin = self.rwin[:irmax] + if with_chiq: + outq = self.fftr(outr) + qmax_out = self.kmax+2.0 + group.q = np.linspace(0, qmax_out, int(1.05 + qmax_out/self.kstep)) + nq = len(group.q) + group.chiq = outq[:nq] + group.chiq_mag = sqrt(outq.real**2 + outq.imag**2)[:nq] + group.chiq_re = outq.real[:nq] + group.chiq_im = outq.imag[:nq] + group.chiq_pha = complex_phase(outq[:nq]) + def get_kweight(self): "if kweight is a list/tuple, use only the first one here" if isinstance(self.kweight, Iterable): @@ -548,7 +559,7 @@ def _residual(self, params, data_only=False, **kws): else: cwt = trans.cwt(diff/eps_k, kweight=trans.kweight) return realimag(cwt).ravel() - else: # 'r' space + else: # 'r' or 'q' space out = [] if all_kweights: chir = [trans.fftf(diff, kweight=kw) for kw in trans.kweight] @@ -562,7 +573,7 @@ def _residual(self, params, data_only=False, **kws): for i, chir_ in enumerate(chir): chir_ = chir_ / (eps_r[i]) out.append(realimag(chir_[irmin:irmax])) - else: + else: # 'q' space chiq = [trans.fftr(c)/eps for c, eps in zip(chir, eps_r)] iqmin = int(max(0, 0.01 + trans.kmin/trans.kstep)) iqmax = int(min(trans.nfft/2, 0.01 + trans.kmax/trans.kstep)) @@ -570,10 +581,11 @@ def _residual(self, params, data_only=False, **kws): out.append( realimag(chiq_[iqmin:iqmax])[::2]) return np.concatenate(out) - def save_outputs(self, rmax_out=10, path_outputs=True): + def save_outputs(self, rmax_out=10, path_outputs=True, with_chiq=True): "save fft outputs, and may also map a refined _bkg to the data chi(k) arrays" def xft(dgroup): - self.transform._xafsft(dgroup.chi, group=dgroup, rmax_out=rmax_out) + self.transform._xafsft(dgroup.chi, group=dgroup, + rmax_out=rmax_out, with_chiq=with_chiq) xft(self.data) xft(self.model) if self.refine_bkg: diff --git a/larch/xafs/pre_edge.py b/larch/xafs/pre_edge.py index 1d8d03d86..86159fb61 100644 --- a/larch/xafs/pre_edge.py +++ b/larch/xafs/pre_edge.py @@ -122,7 +122,7 @@ def _finde0(energy, mu_input, estep=None, use_smooth=True): def flat_resid(pars, en, mu): return pars['c0'] + en * (pars['c1'] + en * pars['c2']) - mu -def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, +def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, npre=1, pre1=None, pre2=None, norm1=None, norm2=None): """pre edge subtraction, normalization for XAFS (straight python) @@ -141,6 +141,7 @@ def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, step: edge jump. If None, it will be determined here. pre1: low E range (relative to E0) for pre-edge fit pre2: high E range (relative to E0) for pre-edge fit + npre: degree for pre-edge fit: 1 for linear, 0 for constant nvict: energy exponent to use for pre-edg fit. See Note norm1: low E range (relative to E0) for post-edge fit norm2: high E range (relative to E0) for post-edge fit @@ -157,10 +158,12 @@ def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, Notes ----- - 1 pre_edge: a line is fit to mu(energy)*energy**nvict over the region, - energy=[e0+pre1, e0+pre2]. pre1 and pre2 default to None, which will set + 1 pre_edge: if npre=1 (default), a line is fit to mu(energy)*energy**nvict + over the region, energy=[e0+pre1, e0+pre2]. + pre1 and pre2 default to None, which will set pre1 = e0 - 2nd energy point, rounded to 5 eV pre2 = roughly pre1/3.0, rounded to 5 eV + if npre=0, a constont (mean value of mu[e0+pre1, e0+pre2] will be used. 2 post-edge: a polynomial of order nnorm is fit to mu(energy)*energy**nvict between energy=[e0+norm1, e0+norm2]. nnorm, norm1, norm2 default to None, @@ -192,7 +195,7 @@ def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, pre1, pre2 = pre2, pre1 ipre1 = index_of(energy-e0, pre1) ipre2 = index_of(energy-e0, pre2) - if ipre2 < ipre1 + 2 + nvict: + if npre==1 and ipre2 < ipre1 + 2 + nvict: pre2 = (energy-e0)[int(ipre1 + 2 + nvict)] if norm2 is None: @@ -215,15 +218,24 @@ def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, # preedge p1 = index_of(energy, pre1+e0) p2 = index_nearest(energy, pre2+e0) - if p2-p1 < 2: - p2 = min(len(energy), p1 + 2) + if npre == 0: + if p2 == p1: + p2 = p2 + 1 + mu_mean = mu[p1:p2].mean() + pre_edge = mu_mean * np.ones(len(energy)) + precoefs = [mu_mean, 0.0] + else: + if p2-p1 < 2: + p2 = min(len(energy), p1 + 2) + omu = mu*energy**nvict + ex = remove_nans(energy[p1:p2], interp=True) + mx = remove_nans(omu[p1:p2], interp=True) - omu = mu*energy**nvict - ex = remove_nans(energy[p1:p2], interp=True) - mx = remove_nans(omu[p1:p2], interp=True) + precoefs = polyfit(ex, mx, 1) + if len(precoefs) == 1: + precoefs.append(0.0) + pre_edge = (precoefs[0] + energy*precoefs[1]) * energy**(-nvict) - precoefs = polyfit(ex, mx, 1) - pre_edge = (precoefs[0] + energy*precoefs[1]) * energy**(-nvict) # normalization p1 = index_of(energy, norm1+e0) p2 = index_nearest(energy, norm2+e0) @@ -252,7 +264,7 @@ def preedge(energy, mu, e0=None, step=None, nnorm=None, nvict=0, pre1=None, @Make_CallArgs(["energy","mu"]) def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, - nvict=0, pre1=None, pre2=None, norm1=None, norm2=None, + nvict=0, npre=1, pre1=None, pre2=None, norm1=None, norm2=None, make_flat=True, _larch=None): """pre edge subtraction, normalization for XAFS @@ -272,6 +284,7 @@ def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, step: edge jump. If None, it will be determined here. pre1: low E range (relative to E0) for pre-edge fit pre2: high E range (relative to E0) for pre-edge fit + npre: degree for pre-edge fit: 1 for linear, 0 for constant nvict: energy exponent to use for pre-edg fit. See Notes. norm1: low E range (relative to E0) for post-edge fit norm2: high E range (relative to E0) for post-edge fit @@ -298,10 +311,12 @@ def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, ----- 1. Supports `First Argument Group` convention, requiring group members `energy` and `mu`. 2. Support `Set XAFS Group` convention within Larch or if `_larch` is set. - 3. pre_edge: a line is fit to mu(energy)*energy**nvict over the region, - energy=[e0+pre1, e0+pre2]. pre1 and pre2 default to None, which will set + 3. pre_edge: if npre=1 (default), a line is fit to mu(energy)*energy**nvict + over the region, energy=[e0+pre1, e0+pre2]. + pre1 and pre2 default to None, which will set pre1 = e0 - 2nd energy point, rounded to 5 eV pre2 = roughly pre1/3.0, rounded to 5 eV + if npre=0, a constont (mean value of mu[e0+pre1, e0+pre2] will be used. 4. post-edge: a polynomial of order nnorm is fit to mu(energy)*energy**nvict between energy=[e0+norm1, e0+norm2]. nnorm, norm1, norm2 default to None, which will set: @@ -330,8 +345,8 @@ def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, if group is not None and e0 is None: e0 = getattr(group, 'e0', None) pre_dat = preedge(energy, mu, e0=e0, step=step, nnorm=nnorm, - nvict=nvict, pre1=pre1, pre2=pre2, norm1=norm1, - norm2=norm2) + nvict=nvict, npre=npre, pre1=pre1, pre2=pre2, + norm1=norm1, norm2=norm2) group = set_xafsGroup(group, _larch=_larch) e0 = pre_dat['e0'] @@ -398,11 +413,11 @@ def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, group.post_edge = pre_dat['post_edge'] group.pre_edge_details = Group() - for attr in ('pre1', 'pre2', 'norm1', 'norm2', 'nnorm', 'nvict'): + for attr in ('pre1', 'pre2', 'norm1', 'norm2', 'nnorm', 'nvict', 'npre'): setattr(group.pre_edge_details, attr, pre_dat.get(attr, None)) - group.pre_edge_details.pre_slope = pre_dat['precoefs'][1] group.pre_edge_details.pre_offset = pre_dat['precoefs'][0] + group.pre_edge_details.pre_slope = pre_dat['precoefs'][1] for i in range(MAX_NNORM): if hasattr(group, 'norm_c%i' % i): @@ -414,7 +429,7 @@ def pre_edge(energy, mu=None, group=None, e0=None, step=None, nnorm=None, group.atsym = getattr(group, 'atsym', None) group.edge = getattr(group, 'edge', None) - if group.atsym is None or group.edge is None: + if (group.atsym is None or group.edge is None) and group.e0 > 10: _atsym, _edge = guess_edge(group.e0) if group.atsym is None: group.atsym = _atsym if group.edge is None: group.edge = _edge diff --git a/larch/xrd/amcsd.py b/larch/xrd/amcsd.py index 5c5616a3a..a5020a6ed 100755 --- a/larch/xrd/amcsd.py +++ b/larch/xrd/amcsd.py @@ -274,7 +274,7 @@ def __init__(self, ams_id=None, ams_db=None, publication=None, mineral=None, def __repr__(self): if self.ams_id is None or self.formula is None: return '' - return f'' + return f'' def get_mineralname(self): minname = self.mineral.name @@ -310,7 +310,7 @@ def ciftext(self): out.append(';') out.append(f"{self.pub_title:s}") out.append(';') - out.append(f"_database_code_amcsd {self.ams_id:07d}") + out.append(f"_database_code_amcsd {self.ams_id}") if self.compound != '': out.append(f"_chemical_compound_source '{self.compound}'") out.append(f"_chemical_formula_sum '{self.formula}'") @@ -564,7 +564,7 @@ def get_pmg_struct(self): ).get_conventional_standard_structure() except: print(f"{err} could not analyze spacegroup for CIF {self.ams_id}") - + def get_unitcell(self): "unitcell as dict, from PMG structure" self.get_pmg_struct() diff --git a/larch/xrd/cif2feff.py b/larch/xrd/cif2feff.py index 99cc7997e..88b21f2ab 100644 --- a/larch/xrd/cif2feff.py +++ b/larch/xrd/cif2feff.py @@ -1,6 +1,7 @@ import os from random import Random from io import StringIO +import numpy as np from xraydb import atomic_symbol, atomic_number, xray_edge from larch.utils import fix_varname, strict_ascii, gformat @@ -32,7 +33,7 @@ def read_cif_structure(ciftext): Arguments --------- - ciftext (string): text of CIF file or name of the CIF file. + ciftext (string): text of CIF file Returns ------- @@ -40,57 +41,198 @@ def read_cif_structure(ciftext): """ if CifParser is None: raise ValueError("CifParser from pymatgen not available. Try 'pip install pymatgen'.") + + if os.path.exists(ciftext): + ciftext = open(ciftext, 'r').read() try: cifstructs = CifParser(StringIO(ciftext), **PMG_CIF_OPTS) parse_ok = True - file_found = True except: parse_ok = False - file_found = False - if os.path.exists(ciftext): - file_found = True - try: - cifstructs = CifParser(ciftext, **PMG_CIF_OPTS) - parse_ok = True - except: - parse_ok = False try: cstruct = cifstructs.parse_structures()[0] except: - raise ValueError('could not get structure from text of CIF file') + raise ValueError('could not get structure from text of CIF') + return cstruct + +def fcompact(val): + """format a float value, removing extra trailing zeros""" + val = f'{val:.6f}' + while val.endswith('0'): + val = val[:-1] + if val.endswith('.'): + val = val + '0' + return val + + +def site_label(site): + coords = ','.join([fcompact(s) for s in site.frac_coords]) + return f'{site.species_string}[{coords}]' + +class CIF_Cluster(): + """ + CIF structure for generating clusters around a specific crystal site, + as used for XAS calculations - if not parse_ok: - if not file_found: - raise FileNotFoundError(f'file {ciftext:s} not found') + """ + def __init__(self, ciftext=None, filename=None, absorber=None, + absorber_site=1, with_h=False, cluster_size=8.0): + self.filename = filename + self.ciftext = ciftext + self.set_absorber(absorber) + self.absorber_site = absorber_site + self.with_h = with_h + self.cluster_size = cluster_size + self.struct = None + + if isinstance(self.absorber, int): + self.absorber = atomic_symbol(self.absorber) + if isinstance(self.absorber, str): + self.absorber_z = atomic_number(self.absorber) + + if ciftext is None and filename is not None: + self.ciftext = open(filename, 'r').read() + + if self.ciftext is not None: + self.parse_ciftext(self.ciftext) + + def set_absorber(self, absorber=None): + self.absorber_z = None + self.absorber = absorber + if isinstance(self.absorber, int): + self.absorber = atomic_symbol(self.absorber) + if isinstance(self.absorber, str): + self.absorber_z = atomic_number(self.absorber) + + def parse_ciftext(self, ciftext=None, absorber=None): + if absorber is not None: + self.set_absorber(absorber) + if ciftext is not None: + self.ciftext = ciftext + self.struct = read_cif_structure(self.ciftext) + self.get_cif_sites() + + def get_cif_sites(self): + """parse sites of CIF structure to get several components: + + struct.sites: list of all sites as parsed by pymatgen + site_labels: list of site labels + unique_sites: list of (site[0], wyckoff sym) for unique xtal sites + unique_map: mapping of all site_labels to unique_site index + absorber_sites: list of unique sites with absorber + + """ + # get equivalent sites, mapping of all sites to unique sites, + # and list of site indexes with absorber + + self.formula = self.struct.composition.reduced_formula + sga = SpacegroupAnalyzer(self.struct) + self.space_group = sga.get_symmetry_dataset().international + + sym_struct = sga.get_symmetrized_structure() + wyckoff_symbols = sym_struct.wyckoff_symbols + + self.site_labels = [] + for site in self.struct.sites: + self.site_labels.append(site_label(site)) + + self.unique_sites = [] + self.unique_map = {} + self.absorber_sites = [] + absorber = '~'*30 if self.absorber is None else self.absorber + for i, sites in enumerate(sym_struct.equivalent_sites): + self.unique_sites.append((sites[0], len(sites), wyckoff_symbols[i])) + for site in sites: + self.unique_map[site_label(site)] = (i+1) + if absorber in site.species_string: + self.absorber_sites.append(i) + + self.atom_sites = {} + self.atom_site_labels = {} + + for i, dat in enumerate(self.unique_sites): + site = dat[0] + label = site_label(site) + for species in site.species: + elem = species.name + if elem in self.atom_sites: + self.atom_sites[elem].append(i+1) + self.atom_site_labels[elem].append(label) + else: + self.atom_sites[elem] = [i+1] + self.atom_site_labels[elem] = [label] + + all_sites = {} + for xat in self.atom_site_labels.keys(): + all_sites[xat] = {} + for i, label in enumerate(self.atom_site_labels[xat]): + all_sites[xat][label] = self.atom_sites[xat][i] + + self.all_sites = all_sites + + + def build_cluster(self, absorber=None, absorber_site=1, cluster_size=None): + if absorber is not None: + self.set_absorber(absorber) + if cluster_size is None: + cluster_size = self.cluster_size + + if absorber_site not in self.atom_sites[self.absorber]: + raise ValueError(f"invalid site for absorber {absorber}: must be in {self.atom_sites[self.absorber]}") + + csize2 = cluster_size**2 + + site_atoms = {} # map xtal site with list of atoms occupying that site + site_tags = {} + + for i, site in enumerate(self.struct.sites): + label = site_label(site) + s_unique = self.unique_map.get(label, 0) + site_species = [e.symbol for e in site.species] + if len(site_species) > 1: + s_els = [s.symbol for s in site.species.keys()] + + s_wts = [s for s in site.species.values()] + site_atoms[i] = rng.choices(s_els, weights=s_wts, k=1000) + site_tags[i] = f'({site.species_string:s})_{s_unique:d}' + else: + site_atoms[i] = [site_species[0]] * 1000 + site_tags[i] = f'{site.species_string:s}_{s_unique:d}' + + # atom0 = self.struct[a_index] + atom0 = self.unique_sites[absorber_site-1][0] + sphere = self.struct.get_neighbors(atom0, self.cluster_size) + + self.symbols = [self.absorber] + self.coords = [[0, 0, 0]] + site0_species = [e.symbol for e in atom0.species] + if len(site0_species) > 1: + self.tags = [f'({atom0.species_string})_{absorber_site:d}'] else: - raise ValueError('invalid text of CIF file') - return cstruct + self.tags = [f'{atom0.species_string}_{absorber_site:d}'] -def cif_sites(ciftext, absorber=None): - "return list of sites for the structure" - cstruct = read_cif_structure(ciftext) - out = cstruct.sites - if absorber is not None: - abname = absorber.lower() - out = [] - for site in cstruct.sites: - species = site.species_string.lower() - if ',' in species and ':' in species: # multi-occupancy site - for siteocc in species.split(','): - sname, occ = siteocc.split(':') - if sname.strip() == abname: - out.append(site) - elif species == abname: - out.append(site) - if len(out) == 0: - out = cstruct.sites[0] - return out + for i, site_dist in enumerate(sphere): + s_index = site_dist[0].index + site_symbol = site_atoms[s_index].pop() + + coords = site_dist[0].coords - atom0.coords + if (coords[0]**2 + coords[1]**2 + coords[2]**2) < csize2: + self.tags.append(site_tags[s_index]) + self.symbols.append(site_symbol) + self.coords.append(coords) + self.molecule = Molecule(self.symbols, self.coords) + +## + +def cif_cluster(ciftext=None, filename=None, absorber=None): + "return list of sites for the structure" + return CIF_Cluster(ciftext=ciftext, filename=filename, absorber=absorber) def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, - site_index=None, extra_titles=None, with_h=False, - version8=True, rng_seed=None): + extra_titles=None, with_h=False, version8=True, rng_seed=None): + """convert CIF text to Feff8 or Feff6l input file Arguments @@ -101,7 +243,6 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, edge (string or None): edge for calculation (see Note 2) [None] cluster_size (float): size of cluster, in Angstroms [8.0] absorber_site (int): index of site for absorber (see Note 3) [1] - site_index (int or None): index of site for absorber (see Note 4) [None] extra_titles (list of str or None): extra title lines to include [None] with_h (bool): whether to include H atoms [False] version8 (bool): whether to write Feff8l input (see Note 5)[True] @@ -122,26 +263,23 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, can be selected by the order in which they are listed in the sites list. This depends on the details of the CIF structure, which can be found with `cif_sites(ciftext)`, starting counting by 1. - 4. to explicitly state the index of the site in the sites list, use - site_index (starting at 1!) 5. if version8 is False, outputs will be written for Feff6l """ - try: - cstruct = read_cif_structure(ciftext) - except ValueError: - return '# could not read CIF file' - global rng if rng_seed is not None: rng.seed(rng_seed) - sgroup = SpacegroupAnalyzer(cstruct).get_symmetry_dataset() - space_group = sgroup["international"] + cluster = CIF_Cluster(ciftext=ciftext, absorber=absorber) + cluster.build_cluster(absorber_site=absorber_site, cluster_size=cluster_size) + + mol = cluster.molecule + + absorber = cluster.absorber + absorber_z = cluster.absorber_z + + # print(cluster.atom_sites, absorber, absorber_site) - if isinstance(absorber, int): - absorber = atomic_symbol(absorber_z) - absorber_z = atomic_number(absorber) if edge is None: edge = 'K' if absorber_z < 58 else 'L3' @@ -150,7 +288,7 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, edge_comment = f'{absorber:s} {edge:s} edge, around {edge_energy:.0f} eV' unique_pot_atoms = [] - for site in cstruct: + for site in cluster.struct: for elem in site.species.elements: if elem.symbol not in unique_pot_atoms: unique_pot_atoms.append(elem.symbol) @@ -163,68 +301,32 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, atlist = ', '.join(atoms_map.keys()) raise ValueError(f'atomic symbol {absorber:s} not listed in CIF data: ({atlist})') - site_atoms = {} # map xtal site with list of atoms occupying that site - site_tags = {} - absorber_count = 0 - for sindex, site in enumerate(cstruct.sites): - site_species = [e.symbol for e in site.species] - if len(site_species) > 1: - s_els = [s.symbol for s in site.species.keys()] - s_wts = [s for s in site.species.values()] - site_atoms[sindex] = rng.choices(s_els, weights=s_wts, k=1000) - site_tags[sindex] = f'({site.species_string:s})_{1+sindex:d}' - else: - site_atoms[sindex] = [site_species[0]] * 1000 - site_tags[sindex] = f'{site.species_string:s}_{1+sindex:d}' - if absorber in site_species: - absorber_count += 1 - if absorber_count == absorber_site: - absorber_index = sindex - - if site_index is not None: - absorber_index = site_index - 1 - - # print("Got sites ", len(cstruct.sites), len(site_atoms), len(site_tags)) - - center = cstruct[absorber_index].coords - sphere = cstruct.get_neighbors(cstruct[absorber_index], cluster_size) - symbols = [absorber] - coords = [[0, 0, 0]] - tags = [f'{absorber:s}_{1+absorber_index:d}'] - - for i, site_dist in enumerate(sphere): - s_index = site_dist[0].index + out_text = ['*** feff input generated by xraylarch cif2feff using pymatgen ***'] - site_symbol = site_atoms[s_index].pop() - tags.append(site_tags[s_index]) - symbols.append(site_symbol) - coords.append(site_dist[0].coords - center) - cluster = Molecule(symbols, coords) - out_text = ['*** feff input generated by xraylarch cif2feff using pymatgen ***'] + out_text.append(f'TITLE Formula: {cluster.formula:s}') + out_text.append(f'TITLE SpaceGroup: {cluster.space_group:s}') if extra_titles is not None: for etitle in extra_titles[:]: - if not etitle.startswith('TITLE '): - etitle = 'TITLE ' + etitle - out_text.append(etitle) - - out_text.append(f'TITLE Formula: {cstruct.composition.reduced_formula:s}') - out_text.append(f'TITLE SpaceGroup: {space_group:s}') - out_text.append(f'TITLE # sites: {cstruct.num_sites}') - - out_text.append('* crystallographics sites: note that these sites may not be unique!') - out_text.append(f'* using absorber at site {1+absorber_index:d} in the list below') - out_text.append(f'* selected as absorber="{absorber:s}", absorber_site={absorber_site:d}') - out_text.append('* index X Y Z species') - for i, site in enumerate(cstruct): + out_text.append('TITLE ' + etitle) + + out_text.append('* ') + out_text.append( '* crystallographic sites:') + out_text.append(f'* To change the absorber site, re-run using `absorber_site`') + out_text.append(f'* with the corresponding site index (counting from 1)') + out_text.append('* site X Y Z Wyckoff species') + + for i, dat in enumerate(cluster.unique_sites): + site, n, wsym = dat fc = site.frac_coords - species_string = fix_varname(site.species_string.strip()) - marker = ' <- absorber' if (i == absorber_index) else '' - out_text.append(f'* {i+1:3d} {fc[0]:.6f} {fc[1]:.6f} {fc[2]:.6f} {species_string:s} {marker:s}') + species_string = site.species_string.strip() + marker = ' <- absorber' if ((i+1) == absorber_site) else '' + s1 = f'{i+1:3d} {fc[0]:.6f} {fc[1]:.6f} {fc[2]:.6f}' + s2 = f'{wsym:>5s} {species_string:s} {marker:s}' + out_text.append(f'* {s1} {s2}') out_text.extend(['* ', '', '']) - if version8: out_text.append(f'EDGE {edge:s}') out_text.append('S02 1.0') @@ -244,15 +346,15 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, out_text.extend(['', 'EXCHANGE 0', '', '* POLARIZATION 0 0 0', '', - 'POTENTIALS', '* IPOT Z Tag']) + 'POTENTIALS', '* IPOT Z Tag']) # loop to find atoms actually in cluster, in case some atom # (maybe fractional occupation) is not included - at_lines = [(0, cluster[0].x, cluster[0].y, cluster[0].z, 0, absorber, tags[0])] + at_lines = [(0, mol[0].x, mol[0].y, mol[0].z, 0, absorber, cluster.tags[0])] ipot_map = {} next_ipot = 0 - for i, site in enumerate(cluster[1:]): + for i, site in enumerate(mol[1:]): sym = site.species_string if sym == 'H' and not with_h: continue @@ -262,19 +364,18 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, next_ipot += 1 ipot_map[sym] = ipot = next_ipot - dist = cluster.get_distance(0, i+1) - at_lines.append((dist, site.x, site.y, site.z, ipot, sym, tags[i+1])) - + dist = mol.get_distance(0, i+1) + at_lines.append((dist, site.x, site.y, site.z, ipot, sym, cluster.tags[i+1])) ipot, z = 0, absorber_z - out_text.append(f' {ipot:4d} {z:4d} {absorber:s}') + out_text.append(f' {ipot:4d} {z:>4d} {absorber:>3s}') for sym, ipot in ipot_map.items(): z = atomic_number(sym) - out_text.append(f' {ipot:4d} {z:4d} {sym:s}') + out_text.append(f' {ipot:4d} {z:>4d} {sym:>3s}') out_text.append('') out_text.append('ATOMS') - out_text.append(f'* x y z ipot tag distance site_info') + out_text.append(f'* x y z ipot tag distance site_info') acount = 0 for dist, x, y, z, ipot, sym, tag in sorted(at_lines, key=lambda x: x[0]): @@ -282,7 +383,7 @@ def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, if acount > 500: break sym = (sym + ' ')[:2] - out_text.append(f' {x: .5f} {y: .5f} {z: .5f} {ipot:4d} {sym:s} {dist:.5f} * {tag:s}') + out_text.append(f' {x:+.5f} {y:+.5f} {z:+.5f} {ipot:2d} {sym:>3s} {dist:.5f} * {tag:s}') out_text.append('') out_text.append('* END') diff --git a/setup.cfg b/setup.cfg index 4d4881d86..64843a0aa 100644 --- a/setup.cfg +++ b/setup.cfg @@ -66,7 +66,6 @@ install_requires = pymatgen<=2024.7.18; python_version <= "3.9" pymatgen>=2024.9.10; python_version > "3.9" mp_api - pycifrw fabio pyfai tabulate