Small cosmetic changes

implemented in ARPACK for calculating a subset of eigenvalues rather than all of
them or just the largest one.

This is invoked by setting algorithm='arnoldi' in the constructor MSSPEC.

    calc = MSSPEC(spectroscopy='EIG', algorithm='arnoldi')

This version finds NEV eigenvalues with the largest magnitudes. The value of
NEV is set in the code to be 2% of the total number. The algorithm also sets
the number of basis vectors, NCV, to be 2*NEV. Options for allowing the user
to set these values as well as other convergence criteria should be added in
later versions.

Files src/msspec/calculator.py and src/msspec/parameters.py have been modified
to accept algorithm='arnoldi'.

A new directory, src/msspec/spec/fortran/eig/ar, has been created and contains
various Fortran files implementing the Arnoldi iteration. Some files are written
in Fortran90 free-format and use some features of the Fortran 2003 standard.

Compilation requires the ARPACK library to be installed and accessible using -larpack

Compilation rules and options have been introduced/modified where necessary in the
files src/options.mk, src/msspec/spec/fortran/Makefile and
src/msspec/spec/fortran/eig_ar.mk

This version has been tested for the CaO example.

Dockerfile

@@ -8,15 +8,13 @@ ARG folder=/opt/msspec user=msspec
 
 # Install system dependencies
 # tools
-RUN apk add bash git make gfortran python3 py3-numpy-f2py 
+RUN apk add bash git make gfortran python3 
 # headers
-RUN apk add python3-dev lapack-dev musl-dev hdf5-dev cairo-dev
+RUN apk add python3-dev lapack-dev musl-dev
 # python packages
-RUN apk add py3-virtualenv py3-pip py3-numpy-dev py3-h5py py3-lxml py3-matplotlib \
+RUN apk add py3-virtualenv py3-pip py3-numpy-dev py3-h5py py3-lxml py3-matplotlib py3-numpy py3-pandas py3-cairo py3-scipy py3-setuptools_scm
-            py3-numpy py3-pandas py3-cairo py3-scipy py3-setuptools_scm \
+RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing py3-wxpython
-            py3-terminaltables ipython
+RUN pip install ase pint terminaltables ipython
-RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/community py3-wxpython
-#RUN pip install ase pint terminaltables ipython
 # for GUI
 RUN apk add ttf-droid adwaita-icon-theme
 RUN apk add build-base
@@ -26,29 +24,19 @@ RUN mkdir -p ${folder}/code
 WORKDIR ${folder}/code
 RUN git clone --branch ${branch} https://${login}:${password}@git.ipr.univ-rennes1.fr/epsi/msspec_python3.git .
 
-RUN virtualenv --system-site-packages ${folder}/.local/src/msspec_venv
-RUN make        pybinding PYTHON=python3 VENV_PATH=${folder}/.local/src/msspec_venv VERBOSE=1
-RUN make -C src sdist     PYTHON=python3 VENV_PATH=${folder}/.local/src/msspec_venv VERBOSE=1
-RUN make -C src frontend  PYTHON=python3 VENV_PATH=${folder}/.local/src/msspec_venv VERBOSE=1
-RUN source ${folder}/.local/src/msspec_venv/bin/activate && pip install src/dist/msspec*tar.gz
-
-
-
 # Build
-#RUN make pybinding NO_VENV=1 PYTHON=python3 VERBOSE=1
+RUN make pybinding NO_VENV=1 PYTHON=python3 VERBOSE=1
-#RUN make -C src sdist PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv 
+RUN make -C src sdist PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv && \
-#&& \
+    pip install src/dist/msspec*tar.gz
-#    pip install src/dist/msspec*tar.gz
 
 # Add a non-privileged user
-#RUN adduser -D -s /bin/bash -h ${folder} ${user}
+RUN adduser -D -s /bin/bash -h ${folder} ${user}
 
 # Set the working directory in the container
-#USER ${user}
+USER ${user}
 
-#RUN virtualenv --system-site-packages ${folder}/.local/src/msspec_venv
+RUN virtualenv --system-site-packages ${folder}/.local/src/msspec_venv
-#RUN source ${folder}/.local/src/msspec_venv/bin/activate && pip install src/dist/msspec*.tar.gz
+RUN make -C src frontend PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv
-#RUN make -C src frontend PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv
 
 
 
@@ -56,7 +44,7 @@ FROM alpine:edge
 # Variables
 ARG folder=/opt/msspec user=msspec
 # Install system dependencies
-RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/community \
+RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing \
 #  hdf5-hl cairo openblas lapack libxml2 libxslt libzlf wxwidgets-gtk3 openjpeg libimagequant \
    nano \
    py3-virtualenv \
@@ -84,8 +72,12 @@ RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/community \
    py3-scipy \
    py3-setuptools_scm \
    py3-wxpython \
-   py3-terminaltables \
+   && pip install \
-   py3-bayesian-optimization \
+   ase \
+   pint \
+   terminaltables \
+   ipython \
+   && pip cache purge \
    # Add a non-privileged user
    && adduser -D -s /bin/bash -h ${folder} ${user}
 
@@ -100,12 +92,10 @@ COPY --from=builder ${folder}/code/src/dist/msspec*tar.gz msspec.tar.gz
 RUN virtualenv --system-site-packages .local/src/msspec_venv && \
     . .local/src/msspec_venv/bin/activate && \
     pip install msspec.tar.gz && \
-    pip install ipython && \
-    pip cache purge && \
     rm -f msspec.tar.gz && \
     mkdir -p .local/bin
 
-COPY --from=builder /root/.local/bin/msspec .local/bin/msspec
+COPY --from=builder ${folder}/.local/bin/msspec .local/bin/msspec
 ENV PATH=${folder}/.local/bin:$PATH
 
 # Run the msspec frontend command on startup

Makefile

@@ -10,8 +10,10 @@ pybinding:
 
 venv:
 ifeq ($(NO_VENV),0)
-	@virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
+#	@virtualenv --python=$(PYTHON_EXE) --prompt="(msspec-$(VERSION)) " $(VENV_PATH)
+	$(PYTHON) -m venv $(VENV_PATH) 
 	$(INSIDE_VENV) python -m ensurepip --upgrade
+	$(INSIDE_VENV) pip install -r src/pip.freeze
 endif
 
 # wget https://bootstrap.pypa.io/get-pip.py && \
@@ -24,7 +26,7 @@ endif
 install: venv pybinding wx
 	@+$(INSIDE_VENV) $(MAKE) -C src sdist
 	@+$(INSIDE_VENV) $(MAKE) -C src frontend
-	@+$(INSIDE_VENV) pip install src/dist/msspec-$(VERSION)*.whl
+	@+$(INSIDE_VENV) pip install --force-reinstall src/dist/msspec-$(VERSION)*.whl
 	@echo "Do not forget to check that $(INSTALL_PREFIX)/bin is set in your \$$PATH"
 
 
@@ -37,10 +39,6 @@ light: VENV_PATH = ./_venv
 light: venv
 	@$(INSIDE_VENV) pip install src/
 
-nogui: VENV_PATH = ./_venv
-nogui: venv pybinding
-	@$(INSIDE_VENV) pip install -e src/
-
 
 _attrdict:
 	# Check if virtualenv python version > 3.3.0

src/msspec/calculator.py

@@ -95,6 +95,7 @@ from msspec.parameters import TMatrixParameters
 from msspec.phagen.fortran.libphagen import main as do_phagen
 from msspec.spec.fortran import _eig_mi
 from msspec.spec.fortran import _eig_pw
+from msspec.spec.fortran import _eig_ar
 from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
 from msspec.spec.fortran import _phd_se_noso_nosp_nosym
 from msspec.spec.fortran import _phd_ce_noso_nosp_nosym
@@ -418,6 +419,8 @@ class _MSCALCULATOR(Calculator):
                     do_spec = _eig_mi.run
                 elif self.global_parameters.algorithm == 'power':
                     do_spec = _eig_pw.run
+                elif self.global_parameters.algorithm == 'arnoldi':
+                    do_spec = _eig_ar.run
                 else:
                     LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
                                  "an allowed combination.".format(self.global_parameters.spectroscopy,
@@ -986,8 +989,8 @@ class _EIG(_MSCALCULATOR):
         _MSCALCULATOR.__init__(self, spectroscopy='EIG', algorithm=algorithm,
                                polarization=polarization, dichroism=dichroism,
                                spinpol=spinpol, folder=folder, txt=txt)
-        if algorithm not in ('inversion', 'power'):
+        if algorithm not in ('inversion', 'power', 'arnoldi'):
-            LOGGER.error("Only the 'inversion' or the 'power' algorithms "
+            LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms "
                          "are supported in EIG spectroscopy mode")
             exit(1)
         self.iodata = iodata.Data('EIG Simulation')

src/msspec/iodata.py

@@ -17,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/iodata.py
-# Last modified: Wed, 26 Feb 2025 11:10:17 +0100
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 
 """
@@ -79,26 +79,19 @@ import ase.io
 from ase.io.extxyz import read_xyz, write_xyz
 import h5py
 import numpy as np
+import wx.grid
 from lxml import etree
+from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
 #from matplotlib.backends.backend_wxagg import FigureCanvasWx as FigureCanvas
 from matplotlib.backends.backend_agg import FigureCanvasAgg
-#from matplotlib.backends.backend_cairo import FigureCanvasCairo as FigureCanvasAgg
+from matplotlib.backends.backend_wxagg import NavigationToolbar2WxAgg
 from matplotlib.figure import Figure
 from terminaltables import AsciiTable
 
 import msspec
+from msspec.msspecgui.msspec.gui.clusterviewer import ClusterViewer
 from msspec.misc import LOGGER
 
-try:
-    import wx.grid
-    from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
-    from matplotlib.backends.backend_wxagg import NavigationToolbar2WxAgg
-    from msspec.msspecgui.msspec.gui.clusterviewer import ClusterViewer
-    has_gui = True
-except ImportError:
-    LOGGER.warning('No modules for GUI')
-    has_gui = False
-
 
 def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1):
     # mix the values of existing theta and new theta and return the
@@ -803,17 +796,11 @@ class Data(object):
         """Pops up a grphical window to show all the defined views of the Data object.
 
         """
-        if has_gui:
+        app = wx.App(False)
-            app = wx.App(False)
+        app.SetAppName('MsSpec Data Viewer')
-            app.SetAppName('MsSpec Data Viewer')
+        frame = _DataWindow(self)
-            frame = _DataWindow(self)
+        frame.Show(True)
-            frame.Show(True)
+        app.MainLoop()
-            app.MainLoop()
-        else:
-            print('**** INFORMATION ****')
-            print('You can not use the Data.view() method since ther is no')
-            print('graphical user interface available in this version of MsSpec.')
-            print("Install WxPython if you need it or use Data.export(...) method instead.")
 
 
 class _DataSetView(object):
@@ -898,18 +885,15 @@ class _DataSetView(object):
                         R = np.sin(np.radians(theta))
                         R_ticks = np.sin(np.radians(theta_ticks))
                     elif proj == 'stereo':
-                        #R = 2 * np.tan(np.radians(theta/2.))
+                        R = 2 * np.tan(np.radians(theta/2.))
-                        #R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
+                        R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
-                        R = theta/90.
-                        R_ticks = theta_ticks/90.
                     #R = np.tan(np.radians(theta/2.))
                     X, Y = np.meshgrid(np.radians(phi), R)
-                    im = axes.pcolormesh(X, Y, Xsec, shading='gouraud')
+                    im = axes.pcolormesh(X, Y, Xsec)
                     axes.set_yticks(R_ticks)
                     axes.set_yticklabels(theta_ticks)
 
-                    cbar = figure.colorbar(im)
+                    figure.colorbar(im)
-                    #im.set_clim(0, 0.0275)
 
                 elif proj == 'polar':
                     values[0] = np.radians(values[0])
@@ -932,7 +916,6 @@ class _DataSetView(object):
             axes.set_ylabel(opts['ylabel'])
             axes.set_xlim(*opts['xlim'])
             axes.set_ylim(*opts['ylim'])
-            #axes.set_axis_off()
             #axes.set_pickradius(5)
             if label:
                 axes.legend()
@@ -1023,426 +1006,425 @@ class _DataSetView(object):
         s += '\tconditions : %s\n' % str(self._selection_conditions)
         return s
 
-if has_gui:
+class _GridWindow(wx.Frame):
-    class _GridWindow(wx.Frame):
+    def __init__(self, dset, parent=None):
-        def __init__(self, dset, parent=None):
+        title = 'Data: ' + dset.title
-            title = 'Data: ' + dset.title
+        wx.Frame.__init__(self, parent, title=title, size=(640, 480))
-            wx.Frame.__init__(self, parent, title=title, size=(640, 480))
+        self.create_grid(dset)
-            self.create_grid(dset)
 
-        def create_grid(self, dset):
+    def create_grid(self, dset):
-            grid = wx.grid.Grid(self, -1)
+        grid = wx.grid.Grid(self, -1)
-            grid.CreateGrid(len(dset), len(dset.columns()))
+        grid.CreateGrid(len(dset), len(dset.columns()))
-            for ic, c in enumerate(dset.columns()):
+        for ic, c in enumerate(dset.columns()):
-                grid.SetColLabelValue(ic, c)
+            grid.SetColLabelValue(ic, c)
-                for iv, v in enumerate(dset[c]):
+            for iv, v in enumerate(dset[c]):
-                    grid.SetCellValue(iv, ic, str(v))
+                grid.SetCellValue(iv, ic, str(v))
 
-    class _ParametersWindow(wx.Frame):
+class _ParametersWindow(wx.Frame):
-        def __init__(self, dset, parent=None):
+    def __init__(self, dset, parent=None):
-            title = 'Parameters: ' + dset.title
+        title = 'Parameters: ' + dset.title
-            wx.Frame.__init__(self, parent, title=title, size=(400, 480))
+        wx.Frame.__init__(self, parent, title=title, size=(400, 480))
-            self.create_tree(dset)
+        self.create_tree(dset)
 
-        def create_tree(self, dset):
+    def create_tree(self, dset):
-            datatree = {}
+        datatree = {}
-            for p in dset.parameters():
+        for p in dset.parameters():
-                is_hidden = p.get('hidden', "False")
+            is_hidden = p.get('hidden', "False")
-                if is_hidden == "True":
+            if is_hidden == "True":
-                    continue
+                continue
-                group = datatree.get(p['group'], [])
+            group = datatree.get(p['group'], [])
-                #strval = str(p['value'] * p['unit'] if p['unit'] else p['value'])
+            #strval = str(p['value'] * p['unit'] if p['unit'] else p['value'])
-                #group.append("{:s} = {:s}".format(p['name'], strval))
+            #group.append("{:s} = {:s}".format(p['name'], strval))
-                group.append("{} = {} {}".format(p['name'], p['value'], p['unit']))
+            group.append("{} = {} {}".format(p['name'], p['value'], p['unit']))
-                datatree[p['group']] = group
+            datatree[p['group']] = group
 
-            tree = wx.TreeCtrl(self, -1)
+        tree = wx.TreeCtrl(self, -1)
-            root = tree.AddRoot('Parameters')
+        root = tree.AddRoot('Parameters')
 
-            for key in list(datatree.keys()):
+        for key in list(datatree.keys()):
-                item0 = tree.AppendItem(root, key)
+            item0 = tree.AppendItem(root, key)
-                for item in datatree[key]:
+            for item in datatree[key]:
-                    tree.AppendItem(item0, item)
+                tree.AppendItem(item0, item)
-            tree.ExpandAll()
+        tree.ExpandAll()
-            tree.SelectItem(root)
+        tree.SelectItem(root)
 
-    class _DataWindow(wx.Frame):
+class _DataWindow(wx.Frame):
-        def __init__(self, data):
+    def __init__(self, data):
-            assert isinstance(data, (Data, DataSet))
+        assert isinstance(data, (Data, DataSet))
 
-            if isinstance(data, DataSet):
+        if isinstance(data, DataSet):
-                dset = data
+            dset = data
-                data = Data()
+            data = Data()
-                data.first = dset
+            data.first = dset
-            self.data = data
+        self.data = data
-            self._filename = None
+        self._filename = None
-            self._current_dset = None
+        self._current_dset = None
 
-            wx.Frame.__init__(self, None, title="", size=(640, 480))
+        wx.Frame.__init__(self, None, title="", size=(640, 480))
 
-            self.Bind(wx.EVT_CLOSE, self.on_close)
+        self.Bind(wx.EVT_CLOSE, self.on_close)
 
-            # Populate the menu bar
+        # Populate the menu bar
-            self.create_menu()
+        self.create_menu()
 
-            # Create the status bar
+        # Create the status bar
-            statusbar = wx.StatusBar(self, -1)
+        statusbar = wx.StatusBar(self, -1)
-            statusbar.SetFieldsCount(3)
+        statusbar.SetFieldsCount(3)
-            statusbar.SetStatusWidths([-2, -1, -1])
+        statusbar.SetStatusWidths([-2, -1, -1])
-            self.SetStatusBar(statusbar)
+        self.SetStatusBar(statusbar)
 
-            # Add the notebook to hold all graphs
+        # Add the notebook to hold all graphs
-            self.notebooks = {}
+        self.notebooks = {}
-            sizer = wx.BoxSizer(wx.VERTICAL)
+        sizer = wx.BoxSizer(wx.VERTICAL)
-            #sizer.Add(self.notebook)
+        #sizer.Add(self.notebook)
-            self.SetSizer(sizer)
+        self.SetSizer(sizer)
 
-            self.Bind(wx.EVT_NOTEBOOK_PAGE_CHANGED, self.on_page_changed)
+        self.Bind(wx.EVT_NOTEBOOK_PAGE_CHANGED, self.on_page_changed)
 
+        self.create_notebooks()
+
+        self.update_title()
+
+    def create_notebooks(self):
+        for key in list(self.notebooks.keys()):
+            nb = self.notebooks.pop(key)
+            nb.Destroy()
+
+        for dset in self.data:
+            nb = wx.Notebook(self, -1)
+            self.notebooks[dset.title] = nb
+            #self.GetSizer().Add(nb, 1, wx.ALL|wx.EXPAND)
+            self.GetSizer().Add(nb, proportion=1, flag=wx.ALL|wx.EXPAND)
+            for view in dset.views():
+                self.create_page(nb, view)
+
+        self.create_menu()
+
+        self.show_dataset(self.data[0].title)
+
+
+    def create_menu(self):
+        menubar = wx.MenuBar()
+        menu1 = wx.Menu()
+        menu1.Append(110, "Open\tCtrl+O")
+        menu1.Append(120, "Save\tCtrl+S")
+        menu1.Append(130, "Save as...")
+        menu1.Append(140, "Export\tCtrl+E")
+        menu1.AppendSeparator()
+        menu1.Append(199, "Close\tCtrl+Q")
+
+        menu2 = wx.Menu()
+        for i, dset in enumerate(self.data):
+            menu_id = 201 + i
+            menu2.AppendRadioItem(menu_id, dset.title)
+            self.Bind(wx.EVT_MENU, self.on_menu_dataset, id=menu_id)
+
+        self.Bind(wx.EVT_MENU, self.on_open, id=110)
+        self.Bind(wx.EVT_MENU, self.on_save, id=120)
+        self.Bind(wx.EVT_MENU, self.on_saveas, id=130)
+        self.Bind(wx.EVT_MENU, self.on_export, id=140)
+        self.Bind(wx.EVT_MENU, self.on_close, id=199)
+
+
+        menu3 = wx.Menu()
+        menu3.Append(301, "Data")
+        menu3.Append(302, "Cluster")
+        menu3.Append(303, "Parameters")
+
+        self.Bind(wx.EVT_MENU, self.on_viewdata, id=301)
+        self.Bind(wx.EVT_MENU, self.on_viewcluster, id=302)
+        self.Bind(wx.EVT_MENU, self.on_viewparameters, id=303)
+
+        menubar.Append(menu1, "&File")
+        menubar.Append(menu2, "&Datasets")
+        menubar.Append(menu3, "&View")
+        self.SetMenuBar(menubar)
+
+    def on_open(self, event):
+        if self.data.is_dirty():
+            mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do '
+                                          'you wish to save before opening'
+                                          'another file ?'),
+                                   'Warning: Unsaved data',
+                                   wx.YES_NO | wx.ICON_WARNING)
+            if mbx.ShowModal() == wx.ID_YES:
+                self.on_saveas(wx.Event())
+            mbx.Destroy()
+
+        wildcard = "HDF5 files (*.hdf5)|*.hdf5"
+        dlg = wx.FileDialog(
+            self, message="Open a file...", defaultDir=os.getcwd(),
+            defaultFile="", wildcard=wildcard, style=wx.FD_OPEN
+            )
+
+        if dlg.ShowModal() == wx.ID_OK:
+            path = dlg.GetPath()
+            self._filename = path
+            self.data = Data.load(path)
             self.create_notebooks()
+        dlg.Destroy()
+        self.update_title()
 
-            self.update_title()
+    def on_save(self, event):
-
+        if self._filename:
-        def create_notebooks(self):
-            for key in list(self.notebooks.keys()):
-                nb = self.notebooks.pop(key)
-                nb.Destroy()
-
-            for dset in self.data:
-                nb = wx.Notebook(self, -1)
-                self.notebooks[dset.title] = nb
-                #self.GetSizer().Add(nb, 1, wx.ALL|wx.EXPAND)
-                self.GetSizer().Add(nb, proportion=1, flag=wx.ALL|wx.EXPAND)
-                for view in dset.views():
-                    self.create_page(nb, view)
-
-            self.create_menu()
-
-            self.show_dataset(self.data[0].title)
-
-
-        def create_menu(self):
-            menubar = wx.MenuBar()
-            menu1 = wx.Menu()
-            menu1.Append(110, "Open\tCtrl+O")
-            menu1.Append(120, "Save\tCtrl+S")
-            menu1.Append(130, "Save as...")
-            menu1.Append(140, "Export\tCtrl+E")
-            menu1.AppendSeparator()
-            menu1.Append(199, "Close\tCtrl+Q")
-
-            menu2 = wx.Menu()
-            for i, dset in enumerate(self.data):
-                menu_id = 201 + i
-                menu2.AppendRadioItem(menu_id, dset.title)
-                self.Bind(wx.EVT_MENU, self.on_menu_dataset, id=menu_id)
-
-            self.Bind(wx.EVT_MENU, self.on_open, id=110)
-            self.Bind(wx.EVT_MENU, self.on_save, id=120)
-            self.Bind(wx.EVT_MENU, self.on_saveas, id=130)
-            self.Bind(wx.EVT_MENU, self.on_export, id=140)
-            self.Bind(wx.EVT_MENU, self.on_close, id=199)
-
-
-            menu3 = wx.Menu()
-            menu3.Append(301, "Data")
-            menu3.Append(302, "Cluster")
-            menu3.Append(303, "Parameters")
-
-            self.Bind(wx.EVT_MENU, self.on_viewdata, id=301)
-            self.Bind(wx.EVT_MENU, self.on_viewcluster, id=302)
-            self.Bind(wx.EVT_MENU, self.on_viewparameters, id=303)
-
-            menubar.Append(menu1, "&File")
-            menubar.Append(menu2, "&Datasets")
-            menubar.Append(menu3, "&View")
-            self.SetMenuBar(menubar)
-
-        def on_open(self, event):
             if self.data.is_dirty():
-                mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do '
+                self.data.save(self._filename)
-                                            'you wish to save before opening'
+        else:
-                                            'another file ?'),
+            self.on_saveas(event)
-                                    'Warning: Unsaved data',
-                                    wx.YES_NO | wx.ICON_WARNING)
-                if mbx.ShowModal() == wx.ID_YES:
-                    self.on_saveas(wx.Event())
-                mbx.Destroy()
 
-            wildcard = "HDF5 files (*.hdf5)|*.hdf5"
+    def on_saveas(self, event):
-            dlg = wx.FileDialog(
+        overwrite = True
-                self, message="Open a file...", defaultDir=os.getcwd(),
+        wildcard = "HDF5 files (*.hdf5)|*.hdf5|All files (*.*)|*.*"
-                defaultFile="", wildcard=wildcard, style=wx.FD_OPEN
+        dlg = wx.FileDialog(
-                )
+            self, message="Save file as ...", defaultDir=os.getcwd(),
+            defaultFile='{}.hdf5'.format(self.data.title.replace(' ','_')),
+            wildcard=wildcard, style=wx.FD_SAVE)
+        dlg.SetFilterIndex(0)
 
-            if dlg.ShowModal() == wx.ID_OK:
+        if dlg.ShowModal() == wx.ID_OK:
-                path = dlg.GetPath()
+            path = dlg.GetPath()
-                self._filename = path
+            if os.path.exists(path):
-                self.data = Data.load(path)
+                mbx = wx.MessageDialog(self, ('This file already exists. '
-                self.create_notebooks()
+                                              'Do you wish to overwrite it ?'),
-            dlg.Destroy()
+                                       'Warning: File exists',
-            self.update_title()
+                                       wx.YES_NO | wx.ICON_WARNING)
-
-        def on_save(self, event):
-            if self._filename:
-                if self.data.is_dirty():
-                    self.data.save(self._filename)
-            else:
-                self.on_saveas(event)
-
-        def on_saveas(self, event):
-            overwrite = True
-            wildcard = "HDF5 files (*.hdf5)|*.hdf5|All files (*.*)|*.*"
-            dlg = wx.FileDialog(
-                self, message="Save file as ...", defaultDir=os.getcwd(),
-                defaultFile='{}.hdf5'.format(self.data.title.replace(' ','_')),
-                wildcard=wildcard, style=wx.FD_SAVE)
-            dlg.SetFilterIndex(0)
-
-            if dlg.ShowModal() == wx.ID_OK:
-                path = dlg.GetPath()
-                if os.path.exists(path):
-                    mbx = wx.MessageDialog(self, ('This file already exists. '
-                                                'Do you wish to overwrite it ?'),
-                                        'Warning: File exists',
-                                        wx.YES_NO | wx.ICON_WARNING)
-                    if mbx.ShowModal() == wx.ID_NO:
-                        overwrite = False
-                    mbx.Destroy()
-                if overwrite:
-                    self.data.save(path)
-                    self._filename = path
-            dlg.Destroy()
-            self.update_title()
-
-        def on_export(self, event):
-            overwrite = True
-            dlg = wx.DirDialog(
-                self, message="Export data...", defaultPath=os.getcwd(),
-                style=wx.DD_DEFAULT_STYLE)
-
-            if dlg.ShowModal() == wx.ID_OK:
-                path = dlg.GetPath()
-                if os.listdir(path):
-                    mbx = wx.MessageDialog(self,
-                                        ('This folder is not empty. '
-                                            'Exporting tour data here may '
-                                            'overwrite its content. Do you wish '
-                                            'to continue ?'),
-                                        'Warning: Folder is not empty',
-                                        wx.YES_NO | wx.ICON_WARNING)
-                    if mbx.ShowModal() == wx.ID_NO:
-                        overwrite = False
-                    mbx.Destroy()
-                self.data.export(path, overwrite)
-            dlg.Destroy()
-
-        def on_viewdata(self, event):
-            dset = self.data[self._current_dset]
-            frame = _GridWindow(dset, parent=self)
-            frame.Show()
-
-        def on_viewcluster(self, event):
-            win = wx.Frame(None, size=wx.Size(480, 340))
-            cluster_viewer = ClusterViewer(win, size=wx.Size(480, 340))
-
-            dset = self.data[self._current_dset]
-            #s = StringIO()
-            #s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
-            #_s = dset.get_parameter(group='Cluster', name='cluster')['value']
-            #print(_s)
-            # rewind to the begining of the string
-            #s.seek(0)
-            #atoms = ase.io.read(s, format='xyz')
-            atoms = dset.get_cluster()
-            cluster_viewer.set_atoms(atoms, rescale=True, center=True)
-            cluster_viewer.rotate_atoms(0., 180.)
-            cluster_viewer.rotate_atoms(-45., -45.)
-            #cluster_viewer.show_emitter(True)
-            win.Show()
-
-        def on_viewparameters(self, event):
-            dset = self.data[self._current_dset]
-            frame = _ParametersWindow(dset, parent=self)
-            frame.Show()
-
-        def on_close(self, event):
-            if self.data.is_dirty():
-                mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do you '
-                                            'really want to quit ?'),
-                                    'Warning: Unsaved data',
-                                    wx.YES_NO | wx.ICON_WARNING)
                 if mbx.ShowModal() == wx.ID_NO:
-                    mbx.Destroy()
+                    overwrite = False
-                    return
+                mbx.Destroy()
-            self.Destroy()
+            if overwrite:
+                self.data.save(path)
+                self._filename = path
+        dlg.Destroy()
+        self.update_title()
+
+    def on_export(self, event):
+        overwrite = True
+        dlg = wx.DirDialog(
+            self, message="Export data...", defaultPath=os.getcwd(),
+            style=wx.DD_DEFAULT_STYLE)
+
+        if dlg.ShowModal() == wx.ID_OK:
+            path = dlg.GetPath()
+            if os.listdir(path):
+                mbx = wx.MessageDialog(self,
+                                       ('This folder is not empty. '
+                                        'Exporting tour data here may '
+                                        'overwrite its content. Do you wish '
+                                        'to continue ?'),
+                                       'Warning: Folder is not empty',
+                                       wx.YES_NO | wx.ICON_WARNING)
+                if mbx.ShowModal() == wx.ID_NO:
+                    overwrite = False
+                mbx.Destroy()
+            self.data.export(path, overwrite)
+        dlg.Destroy()
+
+    def on_viewdata(self, event):
+        dset = self.data[self._current_dset]
+        frame = _GridWindow(dset, parent=self)
+        frame.Show()
+
+    def on_viewcluster(self, event):
+        win = wx.Frame(None, size=wx.Size(480, 340))
+        cluster_viewer = ClusterViewer(win, size=wx.Size(480, 340))
+
+        dset = self.data[self._current_dset]
+        #s = StringIO()
+        #s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
+        #_s = dset.get_parameter(group='Cluster', name='cluster')['value']
+        #print(_s)
+        # rewind to the begining of the string
+        #s.seek(0)
+        #atoms = ase.io.read(s, format='xyz')
+        atoms = dset.get_cluster()
+        cluster_viewer.set_atoms(atoms, rescale=True, center=True)
+        cluster_viewer.rotate_atoms(0., 180.)
+        cluster_viewer.rotate_atoms(-45., -45.)
+        #cluster_viewer.show_emitter(True)
+        win.Show()
+
+    def on_viewparameters(self, event):
+        dset = self.data[self._current_dset]
+        frame = _ParametersWindow(dset, parent=self)
+        frame.Show()
+
+    def on_close(self, event):
+        if self.data.is_dirty():
+            mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do you '
+                                          'really want to quit ?'),
+                                   'Warning: Unsaved data',
+                                   wx.YES_NO | wx.ICON_WARNING)
+            if mbx.ShowModal() == wx.ID_NO:
+                mbx.Destroy()
+                return
+        self.Destroy()
 
 
-        def on_menu_dataset(self, event):
+    def on_menu_dataset(self, event):
-            menu_id = event.GetId()
+        menu_id = event.GetId()
-            dset_name = self.GetMenuBar().FindItemById(menu_id).GetItemLabelText()
+        dset_name = self.GetMenuBar().FindItemById(menu_id).GetItemLabelText()
-            self.show_dataset(dset_name)
+        self.show_dataset(dset_name)
 
 
-        def show_dataset(self, name):
+    def show_dataset(self, name):
-            for nb in list(self.notebooks.values()):
+        for nb in list(self.notebooks.values()):
-                nb.Hide()
+            nb.Hide()
-            self.notebooks[name].Show()
+        self.notebooks[name].Show()
-            self.Layout()
+        self.Layout()
-            self.update_statusbar()
+        self.update_statusbar()
-            self._current_dset = name
+        self._current_dset = name
 
-        def create_page(self, nb, view):
+    def create_page(self, nb, view):
-            # Get the matplotlib figure
+        # Get the matplotlib figure
-            figure = view.get_figure()
+        figure = view.get_figure()
 
-            # Create a panel
+        # Create a panel
-            p = wx.Panel(nb, -1)
+        p = wx.Panel(nb, -1)
 
-            # Create a matplotlib canvas for the figure
+        # Create a matplotlib canvas for the figure
-            canvas = FigureCanvas(p, -1, figure)
+        canvas = FigureCanvas(p, -1, figure)
-            sizer = wx.BoxSizer(wx.VERTICAL)
+        sizer = wx.BoxSizer(wx.VERTICAL)
 
-            toolbar = NavigationToolbar2WxAgg(canvas)
+        toolbar = NavigationToolbar2WxAgg(canvas)
-            toolbar.Realize()
+        toolbar.Realize()
 
-            #sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
+        #sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
-            sizer.Add(toolbar, proportion=0, flag=wx.ALL|wx.EXPAND)
+        sizer.Add(toolbar, proportion=0, flag=wx.ALL|wx.EXPAND)
-            toolbar.update()
+        toolbar.update()
 
-            #sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
+        #sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
-            sizer.Add(canvas, proportion=1, flag=wx.ALL|wx.EXPAND)
+        sizer.Add(canvas, proportion=1, flag=wx.ALL|wx.EXPAND)
 
-            p.SetSizer(sizer)
+        p.SetSizer(sizer)
-            p.Fit()
+        p.Fit()
-            p.Show()
+        p.Show()
 
-            # MPL events
+        # MPL events
-            figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
+        figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
-            figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
+        figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
 
-            nb.AddPage(p, view.title)
+        nb.AddPage(p, view.title)
-            canvas.draw()
+        canvas.draw()
 
 
-        def OLDcreate_page(self, nb, view):
+    def OLDcreate_page(self, nb, view):
-            opts = view._plotopts
+        opts = view._plotopts
-            p = wx.Panel(nb, -1)
+        p = wx.Panel(nb, -1)
 
-            figure = Figure()
+        figure = Figure()
 
-            axes = None
+        axes = None
-            proj = opts['projection']
+        proj = opts['projection']
-            scale = opts['scale']
+        scale = opts['scale']
-            if proj == 'rectilinear':
+        if proj == 'rectilinear':
-                axes = figure.add_subplot(111, projection='rectilinear')
+            axes = figure.add_subplot(111, projection='rectilinear')
-            elif proj in ('polar', 'ortho', 'stereo'):
+        elif proj in ('polar', 'ortho', 'stereo'):
-                axes = figure.add_subplot(111, projection='polar')
+            axes = figure.add_subplot(111, projection='polar')
 
-            canvas = FigureCanvas(p, -1, figure)
+        canvas = FigureCanvas(p, -1, figure)
-            sizer = wx.BoxSizer(wx.VERTICAL)
+        sizer = wx.BoxSizer(wx.VERTICAL)
 
-            toolbar = NavigationToolbar2WxAgg(canvas)
+        toolbar = NavigationToolbar2WxAgg(canvas)
-            toolbar.Realize()
+        toolbar.Realize()
 
-            sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
+        sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
-            toolbar.update()
+        toolbar.update()
 
-            sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
+        sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
 
-            p.SetSizer(sizer)
+        p.SetSizer(sizer)
-            p.Fit()
+        p.Fit()
-            p.Show()
+        p.Show()
 
 
-            for values, label in zip(view.get_data(), opts['legend']):
+        for values, label in zip(view.get_data(), opts['legend']):
-                # if we have only one column to plot, select a bar graph
+            # if we have only one column to plot, select a bar graph
-                if np.shape(values)[0] == 1:
+            if np.shape(values)[0] == 1:
-                    xvalues = list(range(len(values[0])))
+                xvalues = list(range(len(values[0])))
-                    axes.bar(xvalues, values[0], label=label,
+                axes.bar(xvalues, values[0], label=label,
-                            picker=5)
+                          picker=5)
-                    axes.set_xticks(xvalues)
+                axes.set_xticks(xvalues)
+            else:
+                if proj in ('ortho', 'stereo'):
+                    theta, phi, Xsec = cols2matrix(*values)
+                    theta_ticks = np.arange(0, 91, 15)
+                    if proj == 'ortho':
+                        R = np.sin(np.radians(theta))
+                        R_ticks = np.sin(np.radians(theta_ticks))
+                    elif proj == 'stereo':
+                        R = 2 * np.tan(np.radians(theta/2.))
+                        R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
+                    #R = np.tan(np.radians(theta/2.))
+                    X, Y = np.meshgrid(np.radians(phi), R)
+                    im = axes.pcolormesh(X, Y, Xsec)
+                    axes.set_yticks(R_ticks)
+                    axes.set_yticklabels(theta_ticks)
+
+                    figure.colorbar(im)
+
+                elif proj == 'polar':
+                    values[0] = np.radians(values[0])
+                    axes.plot(*values, label=label, picker=5,
+                              marker=opts['marker'])
                 else:
-                    if proj in ('ortho', 'stereo'):
+                    if scale == 'semilogx':
-                        theta, phi, Xsec = cols2matrix(*values)
+                        pltcmd = axes.semilogx
-                        theta_ticks = np.arange(0, 91, 15)
+                    elif scale == 'semilogy':
-                        if proj == 'ortho':
+                        pltcmd = axes.semilogy
-                            R = np.sin(np.radians(theta))
+                    elif scale == 'log':
-                            R_ticks = np.sin(np.radians(theta_ticks))
+                        pltcmd = axes.loglog
-                        elif proj == 'stereo':
-                            R = 2 * np.tan(np.radians(theta/2.))
-                            R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
-                        #R = np.tan(np.radians(theta/2.))
-                        X, Y = np.meshgrid(np.radians(phi), R)
-                        im = axes.pcolormesh(X, Y, Xsec)
-                        axes.set_yticks(R_ticks)
-                        axes.set_yticklabels(theta_ticks)
-
-                        figure.colorbar(im)
-
-                    elif proj == 'polar':
-                        values[0] = np.radians(values[0])
-                        axes.plot(*values, label=label, picker=5,
-                                marker=opts['marker'])
                     else:
-                        if scale == 'semilogx':
+                        pltcmd = axes.plot
-                            pltcmd = axes.semilogx
+                    pltcmd(*values, label=label, picker=5,
-                        elif scale == 'semilogy':
+                           marker=opts['marker'])
-                            pltcmd = axes.semilogy
+            axes.grid(opts['grid'])
-                        elif scale == 'log':
+            axes.set_title(opts['title'])
-                            pltcmd = axes.loglog
+            axes.set_xlabel(opts['xlabel'])
-                        else:
+            axes.set_ylabel(opts['ylabel'])
-                            pltcmd = axes.plot
+            axes.set_xlim(*opts['xlim'])
-                        pltcmd(*values, label=label, picker=5,
+            axes.set_ylim(*opts['ylim'])
-                            marker=opts['marker'])
+            if label:
-                axes.grid(opts['grid'])
+                axes.legend()
-                axes.set_title(opts['title'])
+        axes.autoscale(enable=opts['autoscale'])
-                axes.set_xlabel(opts['xlabel'])
-                axes.set_ylabel(opts['ylabel'])
-                axes.set_xlim(*opts['xlim'])
-                axes.set_ylim(*opts['ylim'])
-                if label:
-                    axes.legend()
-            axes.autoscale(enable=opts['autoscale'])
 
 
-            # MPL events
+        # MPL events
-            figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
+        figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
-            figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
+        figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
 
-            nb.AddPage(p, view.title)
+        nb.AddPage(p, view.title)
 
 
-        def update_statusbar(self):
+    def update_statusbar(self):
-            sb = self.GetStatusBar()
+        sb = self.GetStatusBar()
-            menu_id = self.GetMenuBar().FindMenu('Datasets')
+        menu_id = self.GetMenuBar().FindMenu('Datasets')
-            menu = self.GetMenuBar().GetMenu(menu_id)
+        menu = self.GetMenuBar().GetMenu(menu_id)
-            for item in menu.GetMenuItems():
+        for item in menu.GetMenuItems():
-                if item.IsChecked():
+            if item.IsChecked():
-                    sb.SetStatusText("%s" % item.GetItemLabelText(), 1)
+                sb.SetStatusText("%s" % item.GetItemLabelText(), 1)
-                    break
+                break
 
-        def update_title(self):
+    def update_title(self):
-            title = "MsSpec Data Viewer"
+        title = "MsSpec Data Viewer"
-            if self.data.title:
+        if self.data.title:
-                title += ": " + self.data.title
+            title += ": " + self.data.title
-            if self._filename:
+        if self._filename:
-                title += " [" + os.path.basename(self._filename) + "]"
+            title += " [" + os.path.basename(self._filename) + "]"
-            self.SetTitle(title)
+        self.SetTitle(title)
 
-        def on_mpl_motion(self, event):
+    def on_mpl_motion(self, event):
-            sb = self.GetStatusBar()
+        sb = self.GetStatusBar()
-            try:
+        try:
-                txt = "[{:.3f}, {:.3f}]".format(event.xdata, event.ydata)
+            txt = "[{:.3f}, {:.3f}]".format(event.xdata, event.ydata)
-                sb.SetStatusText(txt, 2)
+            sb.SetStatusText(txt, 2)
-            except Exception:
+        except Exception:
-                pass
+            pass
 
-        def on_mpl_pick(self, event):
+    def on_mpl_pick(self, event):
-            print(event.artist)
+        print(event.artist)
 
-        def on_page_changed(self, event):
+    def on_page_changed(self, event):
-            self.update_statusbar()
+        self.update_statusbar()
 
 
 

src/msspec/iodata_gi.py

@@ -235,8 +235,8 @@ class DataSet(object):
                          float: '{:<20.10e}', complex: 's'}
         self._formats = ((np.integer, '{:<20d}'),
                          (np.floating, '{:<20.10e}'),
-                         (complex, '({0.real:<.10e} {0.imag:<.10e}j)'),
+                         (np.complex, '({0.real:<.10e} {0.imag:<.10e}j)'),
-                         (bool, '{:s}'),
+                         (np.bool, '{:s}'),
                          (str, '{:s}'))
 
 
@@ -450,13 +450,9 @@ class DataSet(object):
         :return: The cluster
         :rtype: :py:class:`ase.Atoms`
         """
-        p = self.get_parameter(group='Cluster', name='cluster')['value']
         s = StringIO()
         s.write(self.get_parameter(group='Cluster', name='cluster')['value'])
-        s.seek(0)
+        return ase.io.read(s, format='xyz')
-        #return ase.io.read(s, format='xyz')
-        cluster = list(read_xyz(s))[-1]
-        return cluster
 
 
     def select(self, *args, **kwargs):
@@ -789,13 +785,13 @@ class Data(object):
                 dset = output.add_dset(dset_name)
                 dset.notes = fd['DATA'][dset_name].attrs['notes']
                 for h5dset in fd['DATA'][dset_name]:
-                    dset.add_columns(**{h5dset: fd['DATA'][dset_name][h5dset][...]})
+                    dset.add_columns(**{h5dset: fd['DATA'][dset_name][h5dset].value})
 
             try:
                 vfile = LooseVersion(fd['MsSpec viewer metainfo'].attrs['version'])
                 if vfile > LooseVersion(msspec.__version__):
                     raise NameError('File was saved with a more recent format')
-                xml = fd['MsSpec viewer metainfo']['info'][...].tobytes()
+                xml = fd['MsSpec viewer metainfo']['info'].value.tostring()
                 root = etree.fromstring(xml)
                 for elt0 in root.iter('parameters'):
                     dset_name = elt0.attrib['dataset']
@@ -858,7 +854,7 @@ class Data(object):
         #win.show()
         #Gtk.main()
         app = _Application(self)
-        exit_status = app.run()#sys.argv)
+        exit_status = app.run(sys.argv)
         sys.exit(exit_status)
 
 class _Application(Gtk.Application):
@@ -951,8 +947,7 @@ class _DataSetView(object):
             if np.shape(values)[0] == 1:
                 xvalues = list(range(len(values[0])))
                 axes.bar(xvalues, values[0], label=label,
-                         # picker=5
+                          picker=5)
-                         )
                 axes.set_xticks(xvalues)
             else:
                 if proj in ('ortho', 'stereo'):
@@ -966,7 +961,7 @@ class _DataSetView(object):
                         R_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
                     #R = np.tan(np.radians(theta/2.))
                     X, Y = np.meshgrid(np.radians(phi), R)
-                    im = axes.pcolormesh(X, Y, Xsec, shading='gouraud')
+                    im = axes.pcolormesh(X, Y, Xsec)
                     axes.set_yticks(R_ticks)
                     axes.set_yticklabels(theta_ticks)
 
@@ -974,7 +969,7 @@ class _DataSetView(object):
 
                 elif proj == 'polar':
                     values[0] = np.radians(values[0])
-                    axes.plot(*values, label=label, #picker=5,
+                    axes.plot(*values, label=label, picker=5,
                               marker=opts['marker'])
                 else:
                     if scale == 'semilogx':
@@ -985,7 +980,7 @@ class _DataSetView(object):
                         pltcmd = axes.loglog
                     else:
                         pltcmd = axes.plot
-                    pltcmd(*values, label=label, #picker=5,
+                    pltcmd(*values, label=label, picker=5,
                            marker=opts['marker'])
             axes.grid(opts['grid'])
             axes.set_title(opts['title'])
@@ -993,7 +988,6 @@ class _DataSetView(object):
             axes.set_ylabel(opts['ylabel'])
             axes.set_xlim(*opts['xlim'])
             axes.set_ylim(*opts['ylim'])
-            #axes.set_pickradius(5)
             if label:
                 axes.legend()
         axes.autoscale(enable=opts['autoscale'])
@@ -1248,7 +1242,7 @@ class _DataWindow(Gtk.ApplicationWindow):
     def on_close(self, action, param):
         if self.data.is_dirty():
             dlg = Gtk.Dialog(title="Warning: Unsaved data",
-                             transient_for=self, modal=True)
+                             transient_for=self, flags=Gtk.DialogFlags.MODAL)
             dlg.add_buttons(Gtk.STOCK_YES, Gtk.ResponseType.YES,
                             Gtk.STOCK_NO, Gtk.ResponseType.NO)
             dlg.set_default_size(150, 100)
@@ -1480,14 +1474,9 @@ class OLD_DataWindow(wx.Frame):
         cluster_viewer = ClusterViewer(win, size=wx.Size(480, 340))
 
         dset = self.data[self._current_dset]
-        #s = StringIO()
+        s = StringIO()
-        #s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
+        s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
-        #_s = dset.get_parameter(group='Cluster', name='cluster')['value']
+        atoms = ase.io.read(s, format='xyz')
-        #print(_s)
-        # rewind to the begining of the string
-        #s.seek(0)
-        #atoms = ase.io.read(s, format='xyz')
-        atoms = dset.get_cluster()
         cluster_viewer.set_atoms(atoms, rescale=True, center=True)
         cluster_viewer.rotate_atoms(0., 180.)
         cluster_viewer.rotate_atoms(-45., -45.)
@@ -1688,7 +1677,7 @@ class OLD_DataWindow(wx.Frame):
 
 
 if __name__ == "__main__":
-    if False:
+    if True:
         data = Data('all my data')
         dset = data.add_dset('Dataset 0')
         X = np.arange(0, 20)
@@ -1724,7 +1713,6 @@ if __name__ == "__main__":
         view.select('x', 'y')
 
         data.view()
-        exit()
+    #import sys
-    import sys
+    #data = Data.load(sys.argv[1])
-    data = Data.load(sys.argv[1])
+    #data.view()
-    data.view()

src/msspec/looper.py

@@ -17,8 +17,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/looper.py
-# Last modified: Wed, 26 Feb 2025 11:15:54 +0100
+# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
-# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
+# Committed by : sylvain tricot <sylvain.tricot@univ-rennes1.fr>
 
 from collections import OrderedDict
 from functools import partial
@@ -92,8 +92,9 @@ class Sweep:
 
 
 class SweepRange:
-    def __init__(self, *sweeps):
+    def __init__(self, *sweeps, passindex=False):
         self.sweeps = sweeps
+        self.passindex = passindex
         self.index = 0
 
         # First check that sweeps that are linked to another on are all included
@@ -157,15 +158,17 @@ class SweepRange:
                     for s in [sweep,] + children:
                         key, value = s[idx]
                         row[key] = value
-                row['sweep_index'] = i
+                if self.passindex:
+                    row['sweep_index'] = i
             return row
         else:
             raise StopIteration
 
     @property
     def columns(self):
-        cols  = ['sweep_index']
+        cols = [sweep.key for sweep in self.sweeps]
-        cols += [sweep.key for sweep in self.sweeps]
+        if self.passindex:
+            cols.append('sweep_index')
         return cols
 
     @property
@@ -199,27 +202,31 @@ class Looper:
         logger.debug("Pipeline called with {}".format(x))
         return self.pipeline(**x)
 
-    def run(self, *sweeps, ncpu=1, **kwargs):
+    def run(self, *sweeps, ncpu=1, passindex=False):
         logger.info("Loop starts...")
         # prepare the list of inputs
-        sr = SweepRange(*sweeps)
+        sr = SweepRange(*sweeps, passindex=passindex)
         items = sr.items
 
         data = []
 
-        t0 = time.time()
-
         if ncpu == 1:
             # serial processing...
             logger.info("serial processing...")
+            t0 = time.time()
+
             for i, values in enumerate(items):
-                values.update(kwargs)
                 result = self._wrapper(values)
                 data.append(result)
+
+            t1 = time.time()
+            dt = t1 - t0
+            logger.info("Processed {:d} sets of inputs in {:.3f} seconds".format(
+                len(sr), dt))
+
         else:
             # Parallel processing...
             chunksize = 1 #int(nsets/ncpu)
-            [values.update(kwargs) for values in items]
             logger.info(("Parallel processing over {:d} cpu (chunksize={:d})..."
                          "").format(ncpu, chunksize))
             t0 = time.time()
@@ -229,23 +236,21 @@ class Looper:
             pool.close()
             pool.join()
 
-        t1 = time.time()
+            t1 = time.time()
-        dt = t1 - t0
+            dt = t1 - t0
-        logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
+            logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
-                        "").format(len(sr), dt))
+                         "").format(len(sr), dt))
 
         # Create the DataFrame
         dfdata = []
-        columns = sr.columns + list(kwargs.keys()) + ['output',]
+        columns = sr.columns + ['output',]
 
         for i in range(len(sr)):
             row = list(items[i].values())
             row.append(data[i])
             dfdata.append(row)
 
-
         df = pd.DataFrame(dfdata, columns=columns)
-        df = df.drop(columns=['sweep_index'])
 
         self.data = df
 
@@ -254,14 +259,14 @@ class Looper:
         # of corresponding dict of parameters {'keyA': [val0,...valn],
         # 'keyB': [val0,...valn], ...}
 
-        # all_xy = []
+        all_xy = []
-        # for irow, row in df.iterrows():
+        for irow, row in df.iterrows():
-            # all_xy.append(row.output[0])
+            all_xy.append(row.output[0])
-            # all_xy.append(row.output[1])
+            all_xy.append(row.output[1])
-        # parameters = df.to_dict()
+        parameters = df.to_dict()
-        # parameters.pop('output')
+        parameters.pop('output')
 
-        return self.data #all_xy, parameters
+        return all_xy, parameters
 
 
 
@@ -271,16 +276,17 @@ class Looper:
 if __name__ == "__main__":
     import numpy as np
     import time
-    import logging
-
-
-    logging.basicConfig(level=logging.DEBUG)
 
+    logger.setLevel(logging.DEBUG)
 
 
     def bar(**kwargs):
-        i = kwargs.get('sweep_index')
+        return 0
-        return np.linspace(0,i,10)
+
+    def post_process(data):
+        x = data.x.unique()
+        y = data.y.unique()
+
 
     theta = Sweep(key='theta', comments="The polar angle",
                   start=-70, stop=70, num=3,
@@ -308,16 +314,7 @@ if __name__ == "__main__":
 
     looper = Looper()
     looper.pipeline = bar
-    other_kws = {'un':1, 'deux':2}
+    data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4,
-    data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4, **other_kws)
+                      passindex=True)
-
-    # Print the dataframe
     print(data)
-
+    #print(data[data.emitter_plane.eq(0)].theta.unique())
-    # Accessing the parameters and ouput values for a given sweep (e.g the last one)
-    print(looper.data.iloc[-1])
-
-    # Post-process the output values. For example here, the output is a 1D-array,
-    # make the sum of sweeps with 'Sr' emitter
-    X = np.array([ x for x in data[data.emitter == 'Sr'].output]).sum(axis=0)
-    print(X)

src/msspec/parameters.py

@@ -770,7 +770,8 @@ class GlobalParameters(BaseParameters):
             Parameter('algorithm', types=str, allowed_values=('expansion',
                                                               'inversion',
                                                               'correlation',
-                                                              'power'),
+                                                              'power',
+                                                              'arnoldi'),
                       default='expansion', doc="""
                       You can choose the algorithm used for the computation of the scattering path operator.
 
@@ -778,6 +779,7 @@ class GlobalParameters(BaseParameters):
                         - '**expansion**', for the Rehr-Albers series expansion
                         - '**correlation**', for the correlation-expansion algorithm
                         - '**power**', for the power method approximation scheme (only for spectroscopy='EIG')
+                        - '**arnoldi**', for computing multiple eigenvalues using Arnoldi iteration (only for spectroscopy='EIG')
 
                       The series expansion algorithm is well suited for high energy since the number of terms
                       required decreases as the energy increases. The exact solution is obtained by the matrix inversion

src/msspec/phagen/fortran/phagen_2.2_dp/msxas3.inc

@@ -28,7 +28,7 @@ c
       integer fl_, rdx_
 c
       parameter ( rdx_     =  1600,
-     $            lmax_    =  80,
+     $            lmax_    =  50,
      $            npss     =  lmax_ + 2,
      $            fl_       =  2*npss + 1,
      $            nef_     =  10,

src/msspec/phagen/fortran/phagen_2.2_dp/phagen_scf_2.2_dp.f

@@ -14625,7 +14625,7 @@ c      check = .true.
 !
       do
 !
-        if (( r_real ( i ) > r_in ) .or. ( i .ge. size(r_real) )) then
+        if ( r_real ( i ) > r_in ) then
 
           exit
 

src/msspec/spec/fortran/Makefile

@@ -1,6 +1,6 @@
-.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw comp_curve clean
+.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve clean
 
-all: phd_se phd_mi phd_ce eig_mi eig_pw comp_curve
+all: phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve
 
 phd_se:
 	@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
@@ -17,6 +17,9 @@ eig_mi:
 eig_pw:
 	@+$(MAKE) -f eig_pw.mk all
 
+eig_ar:
+	@+$(MAKE) -f eig_ar.mk all
+
 comp_curve:
 	@+$(MAKE) -f comp_curve.mk all
 
@@ -26,4 +29,5 @@ clean::
 	@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@
 	@+$(MAKE) -f eig_mi.mk $@
 	@+$(MAKE) -f eig_pw.mk $@
+	@+$(MAKE) -f eig_ar.mk $@
 	@+$(MAKE) -f comp_curve.mk $@

src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90

@@ -0,0 +1,250 @@
+!==============================================================================!
+ module arnoldi_mod
+!==============================================================================!
+    implicit none
+    private
+!
+    integer, parameter :: sp = kind(1.0)
+    integer, parameter :: dp = kind(1.0d0)
+    integer, parameter :: zp = kind((0.0d0,0.0d0))
+!
+    complex(zp), parameter :: one = complex(1.0d0, 0.0d0)
+    complex(zp), parameter :: zero = complex(0.0d0, 0.0d0)
+!
+! Public procedures
+!
+    public :: arnoldi_iteration 
+!
+! Private data
+!
+! ARPACK's debug common block
+!
+    integer :: logfil = 6, ndigit = -3, mcaupd = 1
+    integer :: mgetv0 = 0, msaupd = 0, msaup2 = 0, msaitr = 0, mseigt = 0,  &
+               msapps = 0, msgets = 0, mseupd = 0, mnaupd = 0, mnaup2 = 0,  &
+               mnaitr = 0, mneigh = 0, mnapps = 0, mngets = 0, mneupd = 0,  &
+               mcaup2 = 0, mcaitr = 0, mceigh = 0, mcapps = 0, mcgets = 0,  &
+               mceupd = 0
+!
+    common /debug/ logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt,  &
+                   msapps, msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh,  &
+                   mnapps, mngets, mneupd, mcaupd, mcaup2, mcaitr, mceigh,  &
+                   mcapps, mcgets, mceupd
+!==============================================================================!
+ contains
+!==============================================================================!
+! Public procedures
+!==============================================================================!
+    subroutine arnoldi_iteration (ndim, a, nev, d)
+!
+! Use ARPACK routines to find a few eigenvalues (lambda) and corresponding 
+! eigenvectors (x) for the standard eigenvalue problem:
+!          
+!                        A*x = lambda*x
+! 
+! where A is a general NDIM by NDIM complex matrix
+!
+! This subroutine is based on an ARPACK test program adapted by Logan Boudet
+! as part of his M1 project in Rennes in 2022.
+! 
+    integer,     intent(in)    :: ndim
+    complex(zp), intent(in)    :: a(:,:)
+    integer,     intent(inout) :: nev           ! number of eigenvalues required
+    complex(zp), intent(out)   :: d(:)          ! vector of required eigenvalues
+!
+! Local data
+!
+    character(1), parameter :: bmat = "I"       ! standarg eigenvalue problem
+    character(2), parameter :: which = "LM"     ! find NEV eigenvalues of 
+!                                               ! largest magnitude
+!
+    integer     :: ido, ierr, info, j, lworkl, nconv, ncv 
+    integer     :: iparam(11)
+    integer     :: ipntr(14)
+    real(dp)    :: tol
+    complex(zp) :: sigma
+    logical     :: rvec
+!
+    real(dp),    allocatable :: rd(:,:)
+    real(dp),    allocatable :: rwork(:)
+    complex(zp), allocatable :: ax(:)
+    complex(zp), allocatable :: resid(:)
+    complex(zp), allocatable :: v(:,:)
+    complex(zp), allocatable :: workd(:)
+    complex(zp), allocatable :: workev(:)
+    complex(zp), allocatable :: workl(:)
+    logical,     allocatable :: select(:)
+!
+! External BLAS/LAPACK functions used
+!
+    real(dp), external :: dznrm2, dlapy2 
+!
+!
+    write(6,*)
+    write(6,*) "----------------- BEGIN OF ARNOLDI ITERATION -----------------"
+!
+! NCV is is the largest number of basis vectors that will be used in the 
+! Implicitly Restarted Arnoldi Process. Work per major iteration is 
+! proportional to NDIM*NCV*NCV.
+! 
+! Note: we must have NCV >= NEV + 2, and preferably NCV >= 2*NEV
+!
+!    ncv = max(ceiling(1.125*nev + 15), nev+3)
+    ncv = 2*nev
+!
+    iparam(11) = 0
+    ipntr(14) = 0
+!
+! stopping criteria; machine precision is used if tol <= 0
+!
+    tol = 0.0_dp
+!
+! Algorithm mode
+!
+    iparam(1) = 1           ! exact shift strategy              
+    iparam(3) = 300         ! maximum number of iterations
+    iparam(7) = 1           ! use mode1 of ZNAUPD
+!
+    lworkl = ncv*(3*ncv + 5)
+    allocate(rwork(ncv))
+    allocate(resid(ndim), v(ndim,ncv))
+    allocate(workd(3*ndim), workev(2*ncv), workl(lworkl))
+    allocate(select(ncv))
+!
+! IDO is the reverse communication parameter used to determine action to be taken
+! on return from ZNAUPD. Its initial value must be 0.
+!
+    ido = 0
+!
+! On entry, INFO == 0 instructs ZNAUPD to use a random starting vector.
+! To specify a particular starting vector, set INFO to a non-zero value.
+! The required startng vector should then be supplied in RESID.
+!    
+    info = 0
+!
+! Main reverse communication loop
+!
+    do
+        call znaupd(ido, bmat, ndim, which, nev, tol, resid, ncv, v, ndim,  &
+                    iparam, ipntr, workd, workl, lworkl, rwork, info)
+
+        if (abs(ido) /= 1) exit
+!
+! Matrix-vector multiplication y = A*x
+! Initial vector x is stored starting at  workd(ipntr(1))
+! The result y should be stored in workd(ipntr(2))
+!
+        call matvec(a, ndim, ndim, workd(ipntr(1)), workd(ipntr(2)))
+!
+    end do
+!
+    if (info < 0) then
+        write(6,*)
+        write(6,*) "Error: znaupd returned info = ", info
+        write(6,*) "Check the documentation of znaupd for more information"
+        write(6,*)
+        stop
+    end if
+!
+! No fatal errors, post-process with ZNEUPD to extract computed eigenvalues.
+! Eigenvectors may be also computed by setting RVEC = .TRUE.)
+!
+    if (info == 1) then
+        write(6,*)
+        write(6,*) "Maximum number of iterations reached."
+        write(6,*)
+    else if (info == 3) then
+        write(6,*)
+        write(6,*) "No shifts could be applied during implicit Arnoldi update"
+        write(6,*) "Try increasing NCV."
+        write(6,*)
+    end if      
+!
+    rvec = .false.
+!
+    call zneupd(rvec, 'A', select, d, v, ndim, sigma, workev, bmat, ndim,    &
+                which, nev, tol, resid, ncv, v, ndim, iparam, ipntr, workd,  &
+                workl, lworkl, rwork, ierr)
+!
+    if (ierr /= 0) then
+        write(6,*)
+        write(6,*) "Error: zneupd returned ierr = ", ierr
+        write(6,*) "Check the documentation of zneupd for more information"
+        write(6,*)
+        stop
+    end if
+!
+! Eigenvalues are returned in the one dimensional array D and if RVEC == .TRUE.
+! the corresponding eigenvectors are returned in the first NCONV == IPARAM(5) 
+! columns of the two dimensional array V 
+! 
+    nconv = iparam(5)
+!
+    if (rvec) then
+!
+! Compute the residual norm ||  A*x - lambda*x || for the NCONV accurately 
+! computed eigenvalues and eigenvectors
+!
+        allocate(rd(ncv,3), ax(ndim))
+!
+        do j = 1, nconv
+            call matvec(a, ndim, ndim, v(:,j), ax)
+            call zaxpy(ndim, -d(j), v(:,j), 1, ax, 1)
+            rd(j,1) = real(d(j), dp)
+            rd(j,2) = aimag(d(j))
+            rd(j,3) = dznrm2(ndim, ax, 1) / dlapy2(rd(j,1), rd(j,2))
+        end do
+!
+        call dmout(6, nconv, 3, rd, 2*nev, -6,  &
+                   "Ritz values (Real, Imag) and relative residuals")
+!
+        deallocate(rd, ax)
+!
+    end if
+!
+    write(6,*)
+    write(6,*) "SUMMARY"
+    write(6,*) "======="
+    write(6,*)
+    write(6,*) "Size of the matrix is ", ndim
+    write(6,*) "The number of Ritz values requested is ", nev
+    write(6,*) "The number of Arnoldi vectors generated (NCV) is ", ncv
+    write(6,*) "Portion of the spectrum: ", which
+    write(6,*) "The number of converged Ritz values is ", nconv 
+    write(6,*) "The number of implicit Arnoldi update iterations taken is ", iparam(3)
+    write(6,*) "The number of OP*x is ", iparam(9)
+    write(6,*) "The convergence criterion is ", tol
+    write(6,*) 
+!
+    nev = nconv
+!   
+    write(6,*) "------------------ END OF ARNOLDI ITERATION ------------------"
+!
+    deallocate(rwork)
+    deallocate(resid, v)
+    deallocate(workd, workev, workl)
+    deallocate(select)
+!
+    return
+    end subroutine arnoldi_iteration
+!==============================================================================!
+! Private procedures
+!==============================================================================!
+    subroutine matvec (a, n, lda, x, y)
+!
+! Compute the matrix-vector product a*x, storing the result in y
+!
+    complex(zp), intent(in)  :: a(:,:)
+    integer,     intent(in)  :: n
+    integer,     intent(in)  :: lda
+    complex(zp), intent(in)  :: x(*)
+    complex(zp), intent(out) :: y(*)
+!
+!
+    call zgemv('n', n, n, one, a, lda, x, 1, zero, y, 1)
+!
+    return
+    end subroutine matvec
+!==============================================================================!
+ end module arnoldi_mod
+!==============================================================================!

src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90

@@ -0,0 +1,291 @@
+!==============================================================================!
+ subroutine eig_mat_ar (je, e_kin)
+!==============================================================================!
+! This subroutine stores the G_o T kernel matrix and computes a subset of the
+! the eigenvalues with largest magnitude using Arnoldi iteration methods from 
+! ARPACK
+!
+    use dim_mod,      only: n_gaunt, nl_m
+    use coor_mod,     only: natyp, ncorr, n_prot, sym_at
+    use outfiles_mod, only: outfile2
+    use outunits_mod, only: iuo1, iuo2
+    use trans_mod,    only: lmax, vk, tl
+!
+    use arnoldi_mod,  only: arnoldi_iteration
+!
+    implicit none
+!
+    integer, parameter :: sp = kind(1.0)
+    integer, parameter :: dp = kind(1.0d0)
+    integer, parameter :: cp = kind((0.0,0.0))
+    integer, parameter :: zp = kind((0.0d0,0.0d0))
+!
+! Subroutine arguments
+!
+    integer,  intent(in) :: je
+    real(sp), intent(in) :: e_kin
+!
+! Local data
+!
+    integer,     parameter :: ibess = 3
+    integer,     parameter :: nprint = 10
+    real(dp),    parameter :: pi = acos(-1.0_dp)
+    real(dp),    parameter :: small = 0.0001_dp
+    complex(zp), parameter :: ic = complex(0.0d0,1.0d0)
+    complex(zp), parameter :: zeroc = complex(0.0d0,0.0d0)
+    complex(zp), parameter :: four_pi_i = 4.0_dp*pi*ic
+!
+    integer     :: jatl, jlin, jtyp, jnum, lj, lmj, mj, nbtypj
+    integer     :: katl, klin, ktyp, knum, lk, lmk, mk, nbtypk
+    integer     :: j, jp, l, lin, l_max, l_min, m
+    integer     :: ndim, n_dot, n_eig, nev, nfin, nltwo, npr, n_xmax
+    real(sp)    :: eig, xj, yj, zj, xmax_l
+    real(dp)    :: attkj, xkj, ykj, zkj, rkj, zdkj, krkj
+    complex(zp) :: expkj, sum_l, tlk
+!
+    integer, save :: iout2, iout3
+!
+    real(sp),     allocatable :: w1(:), w2(:)
+    real(dp),     allocatable :: gnt(:)
+    complex(zp),  allocatable :: hl1(:), sm(:,:), ylm(:,:), w(:)
+    character(:), allocatable :: outfile, path
+!
+!
+    if (je == 1) then
+!
+! Name of second output file where eigenvalues will be written
+!
+        n_dot = index(outfile2, '.')
+        outfile = outfile2(1:n_dot)//'egv'
+        path = outfile2(1:n_dot)//'pth'
+        open(newunit=iout2, file=outfile, status='unknown')
+        open(newunit=iout3, file=path, status='unknown')
+!
+    end if
+!
+! Construction of the multiple scattering kernel matrix G_o T.
+! Elements are stored using a linear index LINJ representing (J,LJ)
+!
+! First compute Go T array dimension
+!
+    jlin = 0
+    do jtyp = 1, n_prot
+        nbtypj = natyp(jtyp)
+        lmj = lmax(jtyp,je)
+        do jnum = 1, nbtypj
+            do lj = 0, lmj
+                do mj = -lj, lj
+                    jlin = jlin + 1
+                end do
+            end do
+        end do
+    end do
+!
+    ndim = jlin
+    write(6,*) "GoT matrix SM has dimension ", ndim
+!
+    allocate(sm(ndim,ndim))
+    sm = zeroc
+!
+    nltwo = 2*nl_m
+    allocate(ylm(0:nltwo, -nltwo:nltwo))
+    ylm = zeroc
+!
+    allocate(hl1(0:nltwo))
+    hl1 = zeroc
+!
+    allocate(gnt(0:n_gaunt))
+    gnt = 0.0_dp
+!
+    jlin = 0
+    do jtyp = 1, n_prot
+        nbtypj = natyp(jtyp)
+        lmj = lmax(jtyp,je)
+        do jnum = 1, nbtypj
+            jatl = ncorr(jnum,jtyp)
+            xj = sym_at(1,jatl)
+            yj = sym_at(2,jatl)
+            zj = sym_at(3,jatl)
+            do lj = 0, lmj
+                do mj = -lj, lj
+                    jlin = jlin + 1
+!
+                    klin=0
+                    do ktyp = 1, n_prot
+                        nbtypk = natyp(ktyp)
+                        lmk = lmax(ktyp,je)
+                        do knum = 1, nbtypk
+                            katl = ncorr(knum,ktyp)
+!
+                            if (katl /= jatl) then
+                                xkj = real(sym_at(1,katl) - xj, dp)
+                                ykj = real(sym_at(2,katl) - yj, dp)
+                                zkj = real(sym_at(3,katl) - zj, dp)
+                                rkj = sqrt(xkj*xkj + ykj*ykj + zkj*zkj)
+                                krkj = real(vk(je), dp)*rkj
+                                attkj = exp(-aimag(cmplx(vk(je)))*rkj)
+                                expkj = (xkj + ic*ykj)/rkj
+                                zdkj = zkj/rkj
+                                call sph_har2(2*nl_m, zdkj, expkj, ylm, lmj+lmk)
+                                call besphe2(lmj+lmk+1, ibess, krkj, hl1)
+                            end if
+!
+                            do lk = 0,lmk
+                                l_min = abs(lk-lj)
+                                l_max = lk + lj
+                                tlk = cmplx(tl(lk,1,ktyp,je))
+                                do mk = -lk, lk
+                                    klin = klin + 1
+! 
+                                    sm(klin,jlin) = zeroc
+                                    sum_l = zeroc
+                                    if (katl /= jatl) then
+                                        call gaunt2(lk, mk, lj, mj, gnt)
+                                        do l = l_min, l_max, 2
+                                            m = mj - mk
+                                            if (abs(m) <= l) then
+                                                sum_l = sum_l + (ic**l)*hl1(l)*ylm(l,m)*gnt(l)
+                                            end if
+                                        end do
+                                        sum_l = sum_l*attkj*four_pi_i
+                                    else
+                                        sum_l = zeroc
+                                    end if
+!
+                                    sm(klin,jlin) = tlk*sum_l
+!
+                                end do
+                            end do
+                        end do
+                    end do
+                end do
+            end do
+        end do
+    end do
+!
+    deallocate(ylm, hl1, gnt)
+!
+! Compute subset of eigenvalues of SM using ARPACK
+! 
+! NEV is the number of eigenvalues required, set to 2% of NDIM
+!
+    nev = ceiling(0.02*ndim)
+    allocate(w(nev))
+!
+    call arnoldi_iteration(ndim, sm, nev, w)
+!
+    deallocate(sm)
+!
+! Save results to filestream for easy access from python
+!
+    call save_eigenvalues(w, nev, e_kin)
+!
+! Write results to OUTFILE on unit IOUT2
+!
+    write(iout2,75)
+    write(iout2,110)
+    write(iout2,80) e_kin
+    write(iout2,110)
+    write(iout2,75)
+    write(iout2,105)
+    write(iout2,75)
+!
+    allocate(w1(nev), w2(nev))
+!
+    n_eig = 0
+    xmax_l = 0.0
+    n_xmax = 0
+    do lin = 1, nev
+        eig = real(abs(w(lin)))
+        write(iout2,100) real(w(lin)), aimag(w(lin)), eig
+        if ((eig-xmax_l) > 0.0001) n_xmax = lin
+        xmax_l = max(xmax_l, eig)
+        w1(lin) = eig
+        if (eig > 1.000) then
+            n_eig = n_eig + 1
+        end if
+    end do
+!
+    write(iout2,75)
+    write(iout2,85) xmax_l
+    write(iout2,90) n_xmax
+    write(iout2,95) w(n_xmax)
+    write(iout2,75)
+!
+! Summarize results in main output file
+!
+    call ordre(nev, w1, nfin, w2)
+!
+    write(iuo1,5)
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,15) w2(1)
+    write(iuo1,20) w2(nfin)
+    write(iuo1,10)
+    write(iuo1,10)
+!
+    if (n_eig >= 1) then
+        if (n_eig == 1) then
+            write(iuo1,25) ndim
+        else
+            write(iuo1,30) n_eig, ndim
+        end if
+    end if
+!
+    write(iuo1,65) n_xmax
+    write(iuo1,70) w(n_xmax)
+    write(iuo1,10)
+    write(iout3,100) real(w(n_xmax)), aimag(w(n_xmax))
+!
+    npr = nprint/5
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,35) 5*npr
+    write(iuo1,10)
+    do jp = 0, npr-1
+        j = 5*jp
+        write(iuo1,40) w2(j+1), w2(j+2), w2(j+3), w2(j+4), w2(j+5)
+    enddo
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,45) w2(1)
+    write(iuo2,*) e_kin, w2(1)
+    if (n_eig == 0) then
+        write(iuo1,50)
+    else
+        write(iuo1,55)
+    end if
+    write(iuo1,10)
+    write(iuo1,10)
+    write(iuo1,60)
+!
+    deallocate(w, w1, w2)
+!
+    return
+!
+  5 format(/,11X,'-----------------  EIGENVALUE ANALYSIS ','-----------------')
+ 10 format(11X,'-',54X,'-')
+ 15 format(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-')
+ 20 format(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-')
+ 25 format(11X,'-',6X,'1 EIGENVALUE IS > 1 OF A TOTAL OF ',I8,6X,'-')
+ 30 format(11X,'-',4X,I5,' EIGENVALUES ARE > 1 OF A TOTAL OF ',I8,2X,'-')
+ 35 format(11X,'-',11X,'THE ',I3,' LARGEST EIGENVALUES ARE :',11X,'-')
+ 40 format(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-')
+ 45 format(11X,'-',4X,'SPECTRAL RADIUS OF THE KERNEL MATRIX : ',F8.5,3X,'-')
+ 50 format(11X,'-',5X,'---> THE MULTIPLE SCATTERING SERIES ','CONVERGES',4X,'-')
+ 55 format(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/,11X,'-',   &
+        18X,'SCATTERING SERIES',19X,'-')
+ 60 format(11X,'----------------------------------------','----------------',/)
+ 65 format(11X,'-',5X,'    LABEL OF LARGEST EIGENVALUE   : ',I5,8X,'-')
+ 70 format(11X,'-',5X,'    LARGEST EIGENVALUE  : ','(',F6.3,',',F6.3,')',8X,'-')
+ 75 format('              ')
+ 80 format('    KINETIC ENERGY : ',F7.2,' eV')
+ 85 format('    LARGEST MODULUS OF EIGENVALUE : ',F6.3)
+ 90 format('    LABEL OF LARGEST EIGENVALUE   : ',I5)
+ 95 format('    LARGEST EIGENVALUE            :  (',F6.3,',',F6.3,')')
+100 format(5X,F9.5,2X,F9.5,2X,F9.5)
+105 format(7X,'EIGENVALUES :',3X,'MODULUS :')
+110 format(2X,'-------------------------------')
+!==============================================================================!
+ end subroutine eig_mat_ar
+!==============================================================================!

src/msspec/spec/fortran/eig/ar/eigdif_ar.f

@@ -0,0 +1,103 @@
+C
+C=======================================================================
+C
+      SUBROUTINE EIGDIF_AR
+C
+C   This subroutine computes some of the eigenvalues of the
+C     multiple scattering matrix using Arnoldi iteration as
+C     implemented in ARPACK.
+C
+C                                      Last modified : 21 June 2023
+C
+C      INCLUDE 'spec.inc'
+      USE DIM_MOD
+      USE CONVTYP_MOD
+      USE COOR_MOD, NTCLU => NATCLU, NTP => NATYP
+      USE DEBWAL_MOD
+      USE EIGEN_MOD, NE => NE_EIG, E0 => E0_EIG, EFIN => EFIN_EIG
+      USE OUTFILES_MOD
+      USE OUTUNITS_MOD
+      USE RESEAU_MOD
+      USE TESTS_MOD
+      USE TRANS_MOD
+      USE VALIN_MOD, E1 => E0, PHLUM => PHILUM
+C
+      COMPLEX IC,ONEC
+      COMPLEX TLT(0:NT_M,4,NATM,NE_M)
+C
+C
+      DATA CONV /0.512314/
+C
+      IC=(0.,1.)
+      ONEC=(1.,0.)
+C
+      OPEN(UNIT=IUO2, FILE=OUTFILE2, STATUS='UNKNOWN')
+C
+C   Loop over the energies
+C
+      DO JE=1,NE
+        IF(NE.GT.1) THEN
+          ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
+        ELSEIF(NE.EQ.1) THEN
+          ECIN=E0
+        ENDIF
+        CALL LPM(ECIN,XLPM,*6)
+        XLPM1=XLPM/A
+        IF(IPRINT.GT.0) WRITE(IUO1,56) A,XLPM1
+        IF(ITL.EQ.0) THEN
+          VK(JE)=SQRT(ECIN+VINT)*CONV*A*ONEC
+          VK2(JE)=CABS(VK(JE)*VK(JE))
+        ENDIF
+        GAMMA=1./(2.*XLPM1)
+        IF(IPOTC.EQ.0) THEN
+          VK(JE)=VK(JE)+IC*GAMMA
+        ENDIF
+        IF(I_MFP.EQ.0) THEN
+          VK(JE)=CMPLX(REAL(VK(JE)))
+          VK2(JE)=CABS(VK(JE)*VK(JE))
+        ENDIF
+        IF(I_VIB.EQ.1) THEN
+          IF(IDCM.GE.1) WRITE(IUO1,22)
+          DO JAT=1,N_PROT
+            IF(IDCM.EQ.0) THEN
+              XK2UJ2=VK2(JE)*UJ2(JAT)
+            ELSE
+              XK2UJ2=VK2(JE)*UJ_SQ(JAT)
+              WRITE(IUO1,23) JAT,UJ_SQ(JAT)*A*A
+            ENDIF
+            CALL DWSPH(JAT,JE,XK2UJ2,TLT,I_VIB)
+            DO LAT=0,LMAX(JAT,JE)
+              TL(LAT,1,JAT,JE)=TLT(LAT,1,JAT,JE)
+            ENDDO
+          ENDDO
+       ENDIF
+C
+C     Eigenvalue calculation
+C
+ckmd      IF(I_PWM.EQ.0) THEN
+ckmd        CALL EIG_MAT_MS(JE,ECIN)
+ckmd      ELSE
+ckmd        CALL SPEC_RAD_POWER(JE,ECIN)
+ckmd      ENDIF
+C
+      call eig_mat_ar(je, ecin)
+C
+C
+C     End of the loop on the energy
+C
+      ENDDO
+      GOTO 7
+C
+   6  WRITE(IUO1,55)
+C
+  22  FORMAT(16X,'INTERNAL CALCULATION OF MEAN SQUARE DISPLACEMENTS',/,
+     &25X,' BY DEBYE UNCORRELATED MODEL:',/)
+  23  FORMAT(21X,'ATOM TYPE ',I5,'  MSD = ',F8.6,' ANG**2')
+  55  FORMAT(///,12X,' <<<<<<<<<<  THIS VALUE OF ILPM IS NOT',
+     &'AVAILABLE  >>>>>>>>>>')
+  56  FORMAT(4X,'LATTICE PARAMETER A = ',F6.3,' ANGSTROEMS',4X,
+     *'MEAN FREE PATH = ',F6.3,' * A',//)
+C
+   7  RETURN
+C
+      END

src/msspec/spec/fortran/eig/ar/main.f

@@ -0,0 +1,22 @@
+      SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        USE DIM_MOD
+        IMPLICIT INTEGER (A-Z)
+CF2PY   INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_ 
+CF2PY   INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_
+CF2PY   INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_
+
+        CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, 
+     &  NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, 
+     &  NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, 
+     &  N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
+
+        CALL DO_MAIN()
+        CALL CLOSE_ALL_FILES()
+
+      END SUBROUTINE
+

src/msspec/spec/fortran/eig_ar.mk

@@ -0,0 +1,14 @@
+memalloc_src               := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f 
+cluster_gen_src            := $(wildcard cluster_gen/*.f)
+common_sub_src             := $(wildcard common_sub/*.f)
+renormalization_src        := $(wildcard renormalization/*.f)
+#eig_common_src             := $(wildcard eig/common/*.f)
+eig_common_src             := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f))
+eig_ar_src                 := $(wildcard eig/ar/*.f)
+eig_ar_src_f90             := $(wildcard eig/ar/*.f90)
+
+SRCS   = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(eig_common_src) $(eig_ar_src_f90) $(eig_ar_src)
+MAIN_F = eig/ar/main.f
+SO     = _eig_ar.so
+
+include ../../../options.mk

src/msspec/utils.py

@@ -19,8 +19,8 @@
 # along with this msspec.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Source file  : src/msspec/utils.py
-# Last modified: Wed, 26 Feb 2025 11:15:03 +0100
+# Last modified: Thu, 06 Oct 2022 18:27:24 +0200
-# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
+# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> 1665073644 +0200
 
 
 """
@@ -468,12 +468,8 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
     # the symbol of your emitter
     symbol = cluster[np.where(cluster.get_tags() == emitter_tag)[0][0]].symbol
 
-    if shape.lower() in ('spherical'):
+    assert (diameter != 0 or planes != 0), \
-        assert (diameter != 0 or planes != 0), \
+        "At least one of diameter or planes parameter must be use."
-            "At least one of diameter or planes parameter must be use."
-    elif shape.lower() in ('cylindrical'):
-        assert (diameter != 0 and planes != 0), \
-            "Diameter and planes parameters must be defined for cylindrical shape."
 
     if diameter == 0:
         # calculate the minimal diameter according to the number of planes
@@ -483,7 +479,6 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
 
     # number of repetition in each direction
     rep = int(3*min_diameter/min(a, c))
-    #print("rep = ", rep)
 
     # repeat the cluster
     cluster = cluster.repeat((rep, rep, rep))
@@ -547,7 +542,7 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
         xplan, yplan = get_xypos(cluster, zplan)
         radius = np.sqrt(xplan**2 + yplan**2 + zplan**2)
 
-        if diameter != 0 and shape in ('spherical'):
+        if diameter != 0:
             assert (radius <= diameter/2), ("The number of planes is too high "
                                             "compared to the diameter.")
         radius = max(radius, diameter/2)
@@ -580,90 +575,3 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
     Atoms.translate(cluster, [0, 0, -ze])   # put the emitter in (0,0,0)
 
     return cluster
-
-def shape_cluster(primitive, emitter_tag=0, emitter_plane=0, diameter=0,
-                          planes=0, shape='spherical'):
-
-    """Creates and returns a cluster based on an Atoms object and some
-    parameters.
-
-    :param cluster: the Atoms object used to create the cluster
-    :type  cluster: Atoms object
-    :param emitter_tag: the tag of your emitter
-    :type  emitter_tag: integer
-    :param diameter: the diameter of your cluster in Angströms
-    :type  diameter: float
-    :param planes: the number of planes of your cluster
-    :type  planes: integer
-    :param emitter_plane: the plane where your emitter will be starting by 0
-                          for the first plane
-    :type  emitter_plane: integer
-
-    See :ref:`hemispherical_cluster_faq` for more informations.
-    """
-    # We need the radius of the cluster and the number of planes
-    if shape.lower() in ('ispherical', 'cylindrical'):
-        assert (nplanes != 0 and diameter != 0), "nplanes and diameter cannot be zero for '{}' shape".format(shape)
-    elif shape.lower() in ('spherical'):
-        if diameter <= 0:
-            # find the diameter based on the number of planes
-            assert planes != 0, "planes should be > 0"
-
-
-    n = 3
-    natoms = 0
-    while True:
-        n += 2
-        cluster = primitive.copy()
-        # Repeat the primitive cell
-        cluster = cluster.repeat((n, n, n))
-        center_cluster(cluster)
-
-        # Find the emitter closest to the origin
-        all_tags = cluster.get_tags()
-        are_emitters = all_tags == emitter_tag
-        _ie = np.linalg.norm(cluster[are_emitters].positions, axis=1).argmin()
-        ie = np.nonzero(are_emitters)[0][_ie]
-        # Translate the cluster to this emitter position
-        cluster.translate(-cluster[ie].position)
-        # cut plane at surface and at bottom
-        all_z = np.unique(cluster.positions[:,2])
-        try:
-            zsurf = all_z[all_z >= 0][emitter_plane]
-        except IndexError:
-            # There are not enough planes above the emitter
-            zsurf = all_z.max()
-        try:
-            zbottom = all_z[all_z <= 0][::-1][planes - (emitter_plane+1)]
-        except IndexError:
-            # There are not enough planes below the emitter
-            zbottom = all_z.min()
-        cluster = cut_plane(cluster, z=[zbottom,zsurf])
-        # spherical shape
-        if shape.lower() in ('spherical'):
-            cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,zsurf))
-        if shape.lower() in ('ispherical'):
-            cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,0))
-        elif shape.lower() in ('cylindrical'):
-            cluster = cut_cylinder(cluster, radius=diameter/2)
-        else:
-            raise NameError("Unknown shape")
-        cluster.set_cell(primitive.cell)
-        if len(cluster) <= natoms:
-            break
-        else:
-            natoms = len(cluster)
-
-
-    return cluster
-
-
-
-
-
-
-
-
-
-
-

src/msspec/version.py

@@ -33,18 +33,16 @@ import subprocess
 
 PKGNAME = 'msspec'
 
-thisfile_path = os.path.abspath(__file__)
-thisfile_dir = os.path.dirname(thisfile_path)
-
 try:
     cmd = ["git describe|sed 's/-\([0-9]\+\)-.*/.dev\\1/g'"]
-    result = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, 
+    result = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, shell=True)
-            shell=True, cwd=thisfile_dir)
     __version__ = result.stdout.decode('utf-8').strip()
     if __version__ == "":
         raise
 except Exception as err:
     try:
+        thisfile_path = os.path.abspath(__file__)
+        thisfile_dir = os.path.dirname(thisfile_path)
         versionfile = os.path.join(thisfile_dir, "./VERSION")
         with open(versionfile, "r") as fd:
             __version__ = fd.readline().strip()

src/options.mk

@@ -1,9 +1,9 @@
-PYTHON   = python
+PYTHON   = python3
 PYMAJ    = 3
 PYMIN    = 5
 
 FC       = gfortran
-F2PY     = f2py --f77exec=$(FC) --f90exec=$(FC)
+F2PY     = f2py3 --f77exec=$(FC) --f90exec=$(FC)
 
 NO_VENV  = 0
 DEBUG    = 0
@@ -31,7 +31,7 @@ IFORT_FFLAGS_DBG =
 ################################################################################
 # F2PY CONFIGURATION                                                           #
 ################################################################################
-F2PYFLAGS     = --opt=-O2 -llapack
+F2PYFLAGS     = --opt=-O2 -llapack -larpack 
 F2PYFLAGS_DBG = --debug-capi --debug
 ################################################################################
 
@@ -103,7 +103,7 @@ endif
 
 FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX))
 
-OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS))))
+OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f90,%.o, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS)))))
 
 .PHONY: clean obj all info
 
@@ -141,6 +141,12 @@ $(BUILDDIR)/%.o: %.f
 	$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
 
 
+$(BUILDDIR)/%.o: %.f90
+	@echo "Compiling $@..."
+	mkdir -p $(basename $@)
+	$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
+
+
 $(SO): $(OBJS) $(MAIN_F)
 	@echo "building Python binding $@..."      
 	mkdir -p $(BUILDDIR)

src/pip.freeze

@@ -10,3 +10,5 @@ pycairo
 scipy
 setuptools-scm
 terminaltables
+wheel
+wxPython