epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Compare commits

base: epsi:master
Branches Tags
epsi:master
epsi:devel
epsi:feature/LEED
epsi:feature/arnoldi
epsi:feature/aed
epsi:1.7.post14
epsi:1.7.post13
epsi:1.7.post12
epsi:1.7.post11
epsi:1.7
epsi:1.6
epsi:1.5
epsi:1.4
epsi:1.3
..
compare: epsi:feature/arnoldi
Branches Tags
epsi:devel
epsi:master
epsi:feature/LEED
epsi:feature/arnoldi
epsi:feature/aed
epsi:1.7.post14
epsi:1.7.post13
epsi:1.7.post12
epsi:1.7.post11
epsi:1.7
epsi:1.6
epsi:1.5
epsi:1.4
epsi:1.3

3 Commits
master ... feature/ar

Author SHA1 Message Date
kmdunseath 6da19024ad Small cosmetic changes 2023-06-29 13:51:23 +02:00
kmdunseath deb350a520 For spectrocopy EIG, added the possibility of using Arnoldi iteration techniques 
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.
 2023-06-29 11:11:49 +02:00
kmdunseath f5a2b9779c Modifications to build msspec-python3 
- Makefile: using python -m venv to create the virtual environment
 - src/options.mk: set PYTHON to python3; add filters and rules to treat .f90 files
 - src/pip.freeze: add packages wheel and wxPython
 2023-06-29 10:50:34 +02:00
20 changed files with 2736 additions and 620 deletions
Show Stats Expand all files Collapse all files

50
Dockerfile
View File

@ -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 \
py3-numpy py3-pandas py3-cairo py3-scipy py3-setuptools_scm \
py3-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
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
RUN apk add --no-cache -X http://dl-cdn.alpinelinux.org/alpine/edge/testing 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 -C src sdist PYTHON=python3 NO_VENV=1 VENV_PATH=${folder}/.local/src/msspec_venv
#&& \
# pip install src/dist/msspec*tar.gz
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 && \
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 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 virtualenv --system-site-packages ${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 \
py3-bayesian-optimization \
&& pip install \
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

10
Makefile
View File

@ -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

7
src/msspec/calculator.py
View File

@ -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'):
LOGGER.error("Only the 'inversion' or the 'power' algorithms "
if algorithm not in ('inversion', 'power', 'arnoldi'):
LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms "
"are supported in EIG spectroscopy mode")
exit(1)
self.iodata = iodata.Data('EIG Simulation')

776
src/msspec/iodata.py
View File

@ -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
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# 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.SetAppName('MsSpec Data Viewer')
frame = _DataWindow(self)
frame.Show(True)
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.")
app = wx.App(False)
app.SetAppName('MsSpec Data Viewer')
frame = _DataWindow(self)
frame.Show(True)
app.MainLoop()
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_ticks = 2 * np.tan(np.radians(theta_ticks/2.))
R = theta/90.
R_ticks = theta_ticks/90.
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, shading='gouraud')
im = axes.pcolormesh(X, Y, Xsec)
axes.set_yticks(R_ticks)
axes.set_yticklabels(theta_ticks)
cbar = figure.colorbar(im)
#im.set_clim(0, 0.0275)
figure.colorbar(im)
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):
def __init__(self, dset, parent=None):
title = 'Data: ' + dset.title
wx.Frame.__init__(self, parent, title=title, size=(640, 480))
self.create_grid(dset)
class _GridWindow(wx.Frame):
def __init__(self, dset, parent=None):
title = 'Data: ' + dset.title
wx.Frame.__init__(self, parent, title=title, size=(640, 480))
self.create_grid(dset)
def create_grid(self, dset):
grid = wx.grid.Grid(self, -1)
grid.CreateGrid(len(dset), len(dset.columns()))
for ic, c in enumerate(dset.columns()):
grid.SetColLabelValue(ic, c)
for iv, v in enumerate(dset[c]):
grid.SetCellValue(iv, ic, str(v))
def create_grid(self, dset):
grid = wx.grid.Grid(self, -1)
grid.CreateGrid(len(dset), len(dset.columns()))
for ic, c in enumerate(dset.columns()):
grid.SetColLabelValue(ic, c)
for iv, v in enumerate(dset[c]):
grid.SetCellValue(iv, ic, str(v))
class _ParametersWindow(wx.Frame):
def __init__(self, dset, parent=None):
title = 'Parameters: ' + dset.title
wx.Frame.__init__(self, parent, title=title, size=(400, 480))
self.create_tree(dset)
class _ParametersWindow(wx.Frame):
def __init__(self, dset, parent=None):
title = 'Parameters: ' + dset.title
wx.Frame.__init__(self, parent, title=title, size=(400, 480))
self.create_tree(dset)
def create_tree(self, dset):
datatree = {}
for p in dset.parameters():
is_hidden = p.get('hidden', "False")
if is_hidden == "True":
continue
group = datatree.get(p['group'], [])
#strval = str(p['value'] * p['unit'] if p['unit'] else p['value'])
#group.append("{:s} = {:s}".format(p['name'], strval))
group.append("{} = {} {}".format(p['name'], p['value'], p['unit']))
datatree[p['group']] = group
def create_tree(self, dset):
datatree = {}
for p in dset.parameters():
is_hidden = p.get('hidden', "False")
if is_hidden == "True":
continue
group = datatree.get(p['group'], [])
#strval = str(p['value'] * p['unit'] if p['unit'] else p['value'])
#group.append("{:s} = {:s}".format(p['name'], strval))
group.append("{} = {} {}".format(p['name'], p['value'], p['unit']))
datatree[p['group']] = group
tree = wx.TreeCtrl(self, -1)
root = tree.AddRoot('Parameters')
tree = wx.TreeCtrl(self, -1)
root = tree.AddRoot('Parameters')
for key in list(datatree.keys()):
item0 = tree.AppendItem(root, key)
for item in datatree[key]:
tree.AppendItem(item0, item)
tree.ExpandAll()
tree.SelectItem(root)
for key in list(datatree.keys()):
item0 = tree.AppendItem(root, key)
for item in datatree[key]:
tree.AppendItem(item0, item)
tree.ExpandAll()
tree.SelectItem(root)
class _DataWindow(wx.Frame):
def __init__(self, data):
assert isinstance(data, (Data, DataSet))
class _DataWindow(wx.Frame):
def __init__(self, data):
assert isinstance(data, (Data, DataSet))
if isinstance(data, DataSet):
dset = data
data = Data()
data.first = dset
self.data = data
self._filename = None
self._current_dset = None
if isinstance(data, DataSet):
dset = data
data = Data()
data.first = dset
self.data = data
self._filename = 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
self.create_menu()
# Populate the menu bar
self.create_menu()
# Create the status bar
statusbar = wx.StatusBar(self, -1)
statusbar.SetFieldsCount(3)
statusbar.SetStatusWidths([-2, -1, -1])
self.SetStatusBar(statusbar)
# Create the status bar
statusbar = wx.StatusBar(self, -1)
statusbar.SetFieldsCount(3)
statusbar.SetStatusWidths([-2, -1, -1])
self.SetStatusBar(statusbar)
# Add the notebook to hold all graphs
self.notebooks = {}
sizer = wx.BoxSizer(wx.VERTICAL)
#sizer.Add(self.notebook)
self.SetSizer(sizer)
# Add the notebook to hold all graphs
self.notebooks = {}
sizer = wx.BoxSizer(wx.VERTICAL)
#sizer.Add(self.notebook)
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 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):
def on_save(self, event):
if self._filename:
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()
self.data.save(self._filename)
else:
self.on_saveas(event)
wildcard = "HDF5 files (*.hdf5)|*.hdf5"
dlg = wx.FileDialog(
self, message="Open a file...", defaultDir=os.getcwd(),
defaultFile="", wildcard=wildcard, style=wx.FD_OPEN
)
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()
self._filename = path
self.data = Data.load(path)
self.create_notebooks()
dlg.Destroy()
self.update_title()
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 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:
mbx.Destroy()
return
self.Destroy()
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()
return
self.Destroy()
def on_menu_dataset(self, event):
menu_id = event.GetId()
dset_name = self.GetMenuBar().FindItemById(menu_id).GetItemLabelText()
self.show_dataset(dset_name)
def on_menu_dataset(self, event):
menu_id = event.GetId()
dset_name = self.GetMenuBar().FindItemById(menu_id).GetItemLabelText()
self.show_dataset(dset_name)
def show_dataset(self, name):
for nb in list(self.notebooks.values()):
nb.Hide()
self.notebooks[name].Show()
self.Layout()
self.update_statusbar()
self._current_dset = name
def show_dataset(self, name):
for nb in list(self.notebooks.values()):
nb.Hide()
self.notebooks[name].Show()
self.Layout()
self.update_statusbar()
self._current_dset = name
def create_page(self, nb, view):
# Get the matplotlib figure
figure = view.get_figure()
def create_page(self, nb, view):
# Get the matplotlib figure
figure = view.get_figure()
# Create a panel
p = wx.Panel(nb, -1)
# Create a panel
p = wx.Panel(nb, -1)
# Create a matplotlib canvas for the figure
canvas = FigureCanvas(p, -1, figure)
sizer = wx.BoxSizer(wx.VERTICAL)
# Create a matplotlib canvas for the figure
canvas = FigureCanvas(p, -1, figure)
sizer = wx.BoxSizer(wx.VERTICAL)
toolbar = NavigationToolbar2WxAgg(canvas)
toolbar.Realize()
toolbar = NavigationToolbar2WxAgg(canvas)
toolbar.Realize()
#sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
sizer.Add(toolbar, proportion=0, flag=wx.ALL|wx.EXPAND)
toolbar.update()
#sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
sizer.Add(toolbar, proportion=0, flag=wx.ALL|wx.EXPAND)
toolbar.update()
#sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
sizer.Add(canvas, proportion=1, flag=wx.ALL|wx.EXPAND)
#sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
sizer.Add(canvas, proportion=1, flag=wx.ALL|wx.EXPAND)
p.SetSizer(sizer)
p.Fit()
p.Show()
p.SetSizer(sizer)
p.Fit()
p.Show()
# MPL events
figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
# MPL events
figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
nb.AddPage(p, view.title)
canvas.draw()
nb.AddPage(p, view.title)
canvas.draw()
def OLDcreate_page(self, nb, view):
opts = view._plotopts
p = wx.Panel(nb, -1)
def OLDcreate_page(self, nb, view):
opts = view._plotopts
p = wx.Panel(nb, -1)
figure = Figure()
figure = Figure()
axes = None
proj = opts['projection']
scale = opts['scale']
if proj == 'rectilinear':
axes = figure.add_subplot(111, projection='rectilinear')
elif proj in ('polar', 'ortho', 'stereo'):
axes = figure.add_subplot(111, projection='polar')
axes = None
proj = opts['projection']
scale = opts['scale']
if proj == 'rectilinear':
axes = figure.add_subplot(111, projection='rectilinear')
elif proj in ('polar', 'ortho', 'stereo'):
axes = figure.add_subplot(111, projection='polar')
canvas = FigureCanvas(p, -1, figure)
sizer = wx.BoxSizer(wx.VERTICAL)
canvas = FigureCanvas(p, -1, figure)
sizer = wx.BoxSizer(wx.VERTICAL)
toolbar = NavigationToolbar2WxAgg(canvas)
toolbar.Realize()
toolbar = NavigationToolbar2WxAgg(canvas)
toolbar.Realize()
sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
toolbar.update()
sizer.Add(toolbar, 0, wx.ALL|wx.EXPAND)
toolbar.update()
sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
sizer.Add(canvas, 5, wx.ALL|wx.EXPAND)
p.SetSizer(sizer)
p.Fit()
p.Show()
p.SetSizer(sizer)
p.Fit()
p.Show()
for values, label in zip(view.get_data(), opts['legend']):
# if we have only one column to plot, select a bar graph
if np.shape(values)[0] == 1:
xvalues = list(range(len(values[0])))
axes.bar(xvalues, values[0], label=label,
picker=5)
axes.set_xticks(xvalues)
for values, label in zip(view.get_data(), opts['legend']):
# if we have only one column to plot, select a bar graph
if np.shape(values)[0] == 1:
xvalues = list(range(len(values[0])))
axes.bar(xvalues, values[0], label=label,
picker=5)
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'):
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'])
if scale == 'semilogx':
pltcmd = axes.semilogx
elif scale == 'semilogy':
pltcmd = axes.semilogy
elif scale == 'log':
pltcmd = axes.loglog
else:
if scale == 'semilogx':
pltcmd = axes.semilogx
elif scale == 'semilogy':
pltcmd = axes.semilogy
elif scale == 'log':
pltcmd = axes.loglog
else:
pltcmd = axes.plot
pltcmd(*values, label=label, picker=5,
marker=opts['marker'])
axes.grid(opts['grid'])
axes.set_title(opts['title'])
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'])
pltcmd = axes.plot
pltcmd(*values, label=label, picker=5,
marker=opts['marker'])
axes.grid(opts['grid'])
axes.set_title(opts['title'])
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
figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
# MPL events
figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion)
figure.canvas.mpl_connect('pick_event', self.on_mpl_pick)
nb.AddPage(p, view.title)
nb.AddPage(p, view.title)
def update_statusbar(self):
sb = self.GetStatusBar()
menu_id = self.GetMenuBar().FindMenu('Datasets')
menu = self.GetMenuBar().GetMenu(menu_id)
for item in menu.GetMenuItems():
if item.IsChecked():
sb.SetStatusText("%s" % item.GetItemLabelText(), 1)
break
def update_statusbar(self):
sb = self.GetStatusBar()
menu_id = self.GetMenuBar().FindMenu('Datasets')
menu = self.GetMenuBar().GetMenu(menu_id)
for item in menu.GetMenuItems():
if item.IsChecked():
sb.SetStatusText("%s" % item.GetItemLabelText(), 1)
break
def update_title(self):
title = "MsSpec Data Viewer"
if self.data.title:
title += ": " + self.data.title
if self._filename:
title += " [" + os.path.basename(self._filename) + "]"
self.SetTitle(title)
def update_title(self):
title = "MsSpec Data Viewer"
if self.data.title:
title += ": " + self.data.title
if self._filename:
title += " [" + os.path.basename(self._filename) + "]"
self.SetTitle(title)
def on_mpl_motion(self, event):
sb = self.GetStatusBar()
try:
txt = "[{:.3f}, {:.3f}]".format(event.xdata, event.ydata)
sb.SetStatusText(txt, 2)
except Exception:
pass
def on_mpl_motion(self, event):
sb = self.GetStatusBar()
try:
txt = "[{:.3f}, {:.3f}]".format(event.xdata, event.ydata)
sb.SetStatusText(txt, 2)
except Exception:
pass
def on_mpl_pick(self, event):
print(event.artist)
def on_mpl_pick(self, event):
print(event.artist)
def on_page_changed(self, event):
self.update_statusbar()
def on_page_changed(self, event):
self.update_statusbar()

48
src/msspec/iodata_gi.py
View File

@ -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)'),
(bool, '{:s}'),
(np.complex, '({0.real:<.10e} {0.imag:<.10e}j)'),
(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')
cluster = list(read_xyz(s))[-1]
return cluster
return ase.io.read(s, format='xyz')
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.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()
s = StringIO()
s.write(dset.get_parameter(group='Cluster', name='cluster')['value'])
atoms = ase.io.read(s, format='xyz')
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
data = Data.load(sys.argv[1])
data.view()
#import sys
#data = Data.load(sys.argv[1])
#data.view()

85
src/msspec/looper.py
View File

@ -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
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Mon, 27 Sep 2021 17:49:48 +0200
# 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()
dt = t1 - t0
logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
"").format(len(sr), dt))
t1 = time.time()
dt = t1 - t0
logger.info(("Processed {:d} sets of inputs in {:.3f} seconds"
"").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 = []
# for irow, row in df.iterrows():
# all_xy.append(row.output[0])
# all_xy.append(row.output[1])
# parameters = df.to_dict()
# parameters.pop('output')
all_xy = []
for irow, row in df.iterrows():
all_xy.append(row.output[0])
all_xy.append(row.output[1])
parameters = df.to_dict()
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 np.linspace(0,i,10)
return 0
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, **other_kws)
# Print the dataframe
data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4,
passindex=True)
print(data)
# 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)
#print(data[data.emitter_plane.eq(0)].theta.unique())

4
src/msspec/parameters.py
View File

@ -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

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

@ -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,

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

@ -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

8
src/msspec/spec/fortran/Makefile
View File

@ -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 $@

250
src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90 Normal file
View File

@ -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
!==============================================================================!

1558
src/msspec/spec/fortran/eig/ar/do_main.f Normal file
View File

File diff suppressed because it is too large Load Diff

291
src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 Normal file
View File

@ -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
!==============================================================================!

103
src/msspec/spec/fortran/eig/ar/eigdif_ar.f Normal file
View File

@ -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

22
src/msspec/spec/fortran/eig/ar/main.f Normal file
View File

@ -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

14
src/msspec/spec/fortran/eig_ar.mk Normal file
View File

@ -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

102
src/msspec/utils.py
View File

@ -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
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
# Last modified: Thu, 06 Oct 2022 18:27:24 +0200
# 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), \
"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."
assert (diameter != 0 or planes != 0), \
"At least one of diameter or planes parameter must be use."
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

8
src/msspec/version.py
View File

@ -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,
shell=True, cwd=thisfile_dir)
result = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, shell=True)
__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()

14
src/options.mk
View File

@ -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)

2
src/pip.freeze
View File

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