From cedbfd823e6871350d6ded3fb0357446e274c297 Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Thu, 12 Dec 2019 16:45:12 +0100 Subject: [PATCH] Fix bugs for energy scan in EIG spectroscopy. It was impossible to perform an energy scan in EIG spectroscopy. The calculator and calcio modules were corrected. We can now compute the spectral radius vs energy with and without renormalization for EIG spectroscopy mode The G_n, Sigma_n, Z_n and Lowdin K^2 methods are supported for the inversion algorithm. This commit fixes also several issues in the GUI and in the hemispherical_cluster function --- .gitignore | 4 ++ Makefile | 14 +++--- doc/source/conf.py | 9 +--- doc/source/custom.py | 5 ++- src/Makefile | 6 +++ src/msspec/calcio.py | 9 ++-- src/msspec/calculator.py | 35 ++++++++++++--- src/msspec/iodata.py | 6 +-- src/msspec/parameters.py | 43 ++++++++++++++----- .../spec/fortran/eig/common/eig_mat_ms.f | 2 +- src/msspec/utils.py | 8 ++-- src/msspec/version.py | 2 +- 12 files changed, 98 insertions(+), 45 deletions(-) diff --git a/.gitignore b/.gitignore index e02fe6f..1ced045 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ __pycache__ *.py[cod] *.so +*.o +*.mod *.bak +*.target +htmlcov package/ src/msspec/results.txt diff --git a/Makefile b/Makefile index a7bdeb9..2046bec 100644 --- a/Makefile +++ b/Makefile @@ -15,11 +15,13 @@ endif .PHONY: clean version selfex venv doc -selfex: results purge +selfex: results @echo "Creating the self-extractible setup program... " +# copy the src folder and purge it + @cp -r src src_ + @+$(MAKE) -C ./src_ purge # update the version - @cp ./src/msspec/version.py ./src/msspec/version.py.bak - @echo "__version__ = \"$(VERSION)\"" > ./src/msspec/version.py + @echo "__version__ = \"$(VERSION)\"" > ./src_/msspec/version.py # create the package folder @mkdir -p package # create the *.lsm file @@ -38,10 +40,10 @@ selfex: results purge @echo "Copying-policy: Gnu Library General Public License (GLPL) 2.0" >> msspec.lsm @echo "End" >> msspec.lsm # create the self-extractible archive - @$(MAKESELF) --license "./license.txt" --lsm ./msspec.lsm ./src package/$(SETUPFILE) "Python MsSpec" ./install.sh $(SUPPRESS_OUPUT) -# restore previous version.py file and remove *.lsm file - @mv ./src/msspec/version.py.bak ./src/msspec/version.py + @$(MAKESELF) --license "./license.txt" --lsm ./msspec.lsm ./src_ package/$(SETUPFILE) "Python MsSpec" ./install.sh $(SUPPRESS_OUPUT) +# remove *.lsm file and src_ @rm ./msspec.lsm + @rm -rf ./src_ version: diff --git a/doc/source/conf.py b/doc/source/conf.py index 1c6b515..5f931d1 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -17,13 +17,8 @@ import sys, os # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. -sys.path.insert(0, os.path.abspath('.')) -sys.path.insert(0, os.path.abspath('../src')) -sys.path.insert(0, os.path.abspath('../src/msspec')) -sys.path.insert(0, os.path.abspath('../../src')) -sys.path.insert(0, os.path.abspath('../../src/msspec')) -#sys.path.insert(0, os.path.abspath('../')) - +sys.path.append(os.path.abspath('./')) +sys.path.append(os.path.abspath('../../src')) import custom diff --git a/doc/source/custom.py b/doc/source/custom.py index 1822edb..1c32387 100644 --- a/doc/source/custom.py +++ b/doc/source/custom.py @@ -8,6 +8,10 @@ import re import numpy as np from lxml import etree +print(os.getcwd()) +sys.path.append('../../src') +from msspec.calculator import MSSPEC + def get_package_version(): p = subprocess.run(["git describe|sed 's/-/.post/'|cut -d'-' -f1"], shell=True, stdout=subprocess.PIPE) @@ -85,7 +89,6 @@ def generate_parameters(spectroscopy=None): return content content = "" - from msspec.calculator import MSSPEC c = MSSPEC(spectroscopy='PED' if spectroscopy==None else spectroscopy) c.shutdown() diff --git a/src/Makefile b/src/Makefile index 92b4de1..101317d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -36,12 +36,18 @@ tests: pybinding clean: @echo "Cleaning all..." + @find ./ -type f -name '*.pyc' -exec rm -f {} + + @find ./ -type d -name '__pycache__' -exec rm -rf {} + @+$(MAKE) -C msspec/spec/fortran clean $(SUPPRESS_OUPUT) @+$(MAKE) -C msspec/phagen/fortran clean $(SUPPRESS_OUPUT) # remove previous sdist @rm -rf dist @rm -rf *.egg* +purge: clean + @echo "Removing also shared objects..." + @find ./ -type f -name '*.so' -exec rm -f {} + + help: @echo "help message" diff --git a/src/msspec/calcio.py b/src/msspec/calcio.py index f30ce10..a4c8004 100644 --- a/src/msspec/calcio.py +++ b/src/msspec/calcio.py @@ -589,7 +589,7 @@ class SpecIO(object): line = fillstr(line, str(p.get_parameter('eigval_idamp')), 39, 'left') content += line - if p.get_parameter('calctype_spectro') == "EIG": + if p.calctype_spectro == "EIG": nlines = int(np.ceil(p.eigval_ne / 4.)) else: nlines = 1 @@ -678,10 +678,7 @@ class SpecIO(object): else: ibwd_arr = np.zeros((nat), dtype=np.int) thbwd_arr = np.ones((nat)) - print("nat", nat) - print("thbwd_arr", thbwd_arr) for at in p.extra_atoms: - print(at) i = at.get('proto_index') - 1 thfwd_arr[i] = at.get('forward_angle') thbwd_arr[i] = at.get('backward_angle') @@ -947,8 +944,8 @@ class SpecIO(object): skip = rows2skip[spectro] data = np.loadtxt(filename, skiprows=skip, unpack=True) - if len(data.shape) <= 1: - data = data.reshape((1, data.shape[0])) + #if len(data.shape) <= 1: + # data = data.reshape((1, data.shape[0])) return data def load_facdif(self, filename='facdif1.dat'): diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 325b384..7ba69bd 100644 --- a/src/msspec/calculator.py +++ b/src/msspec/calculator.py @@ -24,6 +24,7 @@ For more information on MsSpec, follow this import os +import shutil import sys import re import inspect @@ -68,6 +69,7 @@ from msspec.spec.fortran import eig_pw from terminaltables.ascii_table import AsciiTable + class _MSCALCULATOR(Calculator): """ This class defines an ASE calculator for doing Multiple scattering @@ -282,15 +284,19 @@ class _MSCALCULATOR(Calculator): self.modified = self.modified or mod1 #or mod0 or mod1 - if self.modified: + #if self.modified: + if mod1: # run phagen #self._make('tmatrix') os.chdir(os.path.join(self.tmp_folder, 'output')) do_phagen() # rename some output files to be more explicit - os.rename('fort.10', 'cluster.clu') - os.rename('fort.35', 'tmatrix.tl') - os.rename('fort.55', 'tmatrix.rad') + #os.rename('fort.10', 'cluster.clu') + #os.rename('fort.35', 'tmatrix.tl') + #os.rename('fort.55', 'tmatrix.rad') + shutil.copy('fort.10', 'cluster.clu') + shutil.copy('fort.35', 'tmatrix.tl') + shutil.copy('fort.55', 'tmatrix.rad') os.chdir(self.init_folder) @@ -487,6 +493,12 @@ class _MSCALCULATOR(Calculator): for p in parameters: _.append(p) return _ + + def add_cluster_to_dset(self, dset): + clusbuf = StringIO() + self.atoms.info['absorber'] = self.atoms.absorber + self.atoms.write(clusbuf, format='xyz') + dset.add_parameter(group='Cluster', name='cluster', value=clusbuf.getvalue(), hidden="True") def shutdown(self): """Removes the temporary folder and all its content. @@ -899,7 +911,7 @@ class _EIG(_MSCALCULATOR): self.get_tmatrix() self.run_spec() - # Now load the data + # Now create a dataset ndset = len(self.iodata) dset = self.iodata.add_dset('Eigen values calculation [{:d}]'.format(ndset)) for p in self.get_parameters(): @@ -911,9 +923,20 @@ class _EIG(_MSCALCULATOR): results_fname = os.path.join(self.tmp_folder, 'output/results.dat') data = self.specio.load_results(results_fname) + + energy, spec_rad = data + dset.add_columns(energy=energy, spectral_radius=spec_rad) - dset.add_columns(energy=data[:,0], eigen_value=data[:,1]) + # create a view + title = 'Spectral radii versus kinetic energy' + xlabel = r'Kinetic energy (eV)' + ylabel = r'Spectral radius' + view = dset.add_view("Spectral radius", title=title, xlabel=xlabel, ylabel=ylabel) + view.select('energy', 'spectral_radius') + # add the cluster to the dataset + self.add_cluster_to_dset(dset) + return self.iodata diff --git a/src/msspec/iodata.py b/src/msspec/iodata.py index b0d8065..393d290 100644 --- a/src/msspec/iodata.py +++ b/src/msspec/iodata.py @@ -429,7 +429,7 @@ class DataSet(object): all_indices = np.arange(0, len(self)) indices = all_indices if len(self) > max_len: - indices = list(range(max_len/2)) + list(range(-max_len/2, 0)) + indices = list(range(int(max_len/2))) + list(range(int(-max_len/2), 0)) _i = 0 for i in indices: @@ -940,7 +940,7 @@ class _DataWindow(wx.Frame): wildcard = "HDF5 files (*.hdf5)|*.hdf5" dlg = wx.FileDialog( self, message="Open a file...", defaultDir=os.getcwd(), - defaultFile="", wildcard=wildcard, style=wx.OPEN + defaultFile="", wildcard=wildcard, style=wx.FD_OPEN ) if dlg.ShowModal() == wx.ID_OK: @@ -964,7 +964,7 @@ class _DataWindow(wx.Frame): dlg = wx.FileDialog( self, message="Save file as ...", defaultDir=os.getcwd(), defaultFile='{}.hdf5'.format(self.data.title.replace(' ','_')), - wildcard=wildcard, style=wx.SAVE) + wildcard=wildcard, style=wx.FD_SAVE) dlg.SetFilterIndex(0) if dlg.ShowModal() == wx.ID_OK: diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 5537539..16f6254 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -1,5 +1,5 @@ # coding: utf-8 -# vim: set et ts=4 sw=4 sts nu fdm=indent: +# vim: set et ts=4 sw=4 nu fdm=indent: """ Module parameters ================= @@ -1383,13 +1383,19 @@ class CalculationParameters(BaseParameters): doc=textwrap.dedent(""" The scattering order. Only meaningful for the 'expansion' algorithm. Its value is limited to 10.""")), - Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n'), + Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n', + 'Z_n', 'K^2'), types=(type(None), str), default=None, doc=textwrap.dedent(""" Enable the calculation of the coefficients for the renormalization of the multiple scattering series. - Only meaningful for the 'expansion' algorithm. You can choose to - renormalize in terms of the Sigma_n or G_n matrices.""")), + You can choose to renormalize in terms of the :math:`\\Sigma_n`, the + :math:`G_n`, the :math:`Z_n` or the Lowdin :math:`K^2` matrices. + + .. note:: + + Please note that the last two renormalization methods only work + for the Eigen Value mode.""")), Parameter('renormalization_omega', types=(int,float,complex), default=1.+0j, doc=textwrap.dedent(""" @@ -1520,14 +1526,29 @@ class CalculationParameters(BaseParameters): LOGGER.info('Scattering order set to %s', p.value) def bind_renormalization_mode(self, p): - if p.value is None: - self.spec_parameters.calc_iren = 0 - else: - if p.value.lower() == 'G_n'.lower(): - self.spec_parameters.calc_iren = 1 - elif p.value.lower() == 'Sigma_n'.lower(): - self.spec_parameters.calc_iren = 2 + iren_map = {None : 0, + 'G_n' : 1, + 'Sigma_n': 2, + 'Z_n' : 3, + 'K^2' : 4} + # Check that the method is neither 'Z_n' nor 'K^2' for other + # 'spetroscopy' than EIG + try: + if (self.global_parameters.spectroscopy == 'PED' + and self.global_parameters.algorithm == 'expansion'): + assert( p.value in (None, 'Sigma_n', 'G_n') ) + elif (self.global_parameters.spectroscopy == 'EIG' + and self.global_parameters.algorithm == 'inversion'): + assert( p.value in p.allowed_values) + else: + assert( p.value is None ) + self.spec_parameters.calc_iren = iren_map[p.value] LOGGER.info("Renormalization activated with \'{}\' method".format(p.value)) + except AssertionError: + LOGGER.error("Renormalization method \'{}\' is not compatible with this " + "spectroscopy/algorithm!".format(p.value)) + exit(1) + def bind_renormalization_omega(self, p): omega = complex(p.value) diff --git a/src/msspec/spec/fortran/eig/common/eig_mat_ms.f b/src/msspec/spec/fortran/eig/common/eig_mat_ms.f index 71b75cc..b29c9ae 100644 --- a/src/msspec/spec/fortran/eig/common/eig_mat_ms.f +++ b/src/msspec/spec/fortran/eig/common/eig_mat_ms.f @@ -249,7 +249,7 @@ CKMD ENDIF CKMD IF (I_REN.NE.0) THEN - WRITE(IUO2,*) E_KIN,W2(1),REN_R,REN_I + WRITE(IUO2,*) E_KIN,W2(1) !,REN_R,REN_I ELSE WRITE(IUO2,*) E_KIN,W2(1) ENDIF diff --git a/src/msspec/utils.py b/src/msspec/utils.py index 82b9367..afd2708 100644 --- a/src/msspec/utils.py +++ b/src/msspec/utils.py @@ -251,12 +251,12 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, p """ def get_xypos(cluster, ze, symbol=None): - nmin = None + nmin = np.inf for atom in cluster: if ze - eps < atom.z < ze + eps and (atom.symbol == symbol or symbol == None): n = np.sqrt(atom.x**2 + atom.y**2) - if (n < nmin) or (nmin is None): + if (n < nmin): nmin = n iatom = atom.index @@ -270,6 +270,8 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, p c = cell[:, 2].max() # a lattice parameter a = cell[:, 0].max() # a lattice parameter p = np.alen(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # the number of planes in the cluster + with np.printoptions(threshold=np.inf): + print(cluster.get_tags() == emitter_tag) symbol = cluster[np.where(cluster.get_tags() == emitter_tag)[0][0]].symbol # the symbol of your emitter assert (diameter != 0 or planes != 0), "At least one of diameter or planes parameter must be use." @@ -329,4 +331,4 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, p ze = all_z[- emitter_plane - 1] # the z-coordinate of the emitter Atoms.translate(cluster, [0, 0, -ze]) # put the emitter in (0,0,0) - return cluster \ No newline at end of file + return cluster diff --git a/src/msspec/version.py b/src/msspec/version.py index bc00e17..cd199df 100644 --- a/src/msspec/version.py +++ b/src/msspec/version.py @@ -1,3 +1,3 @@ # This file is updated at package creation # DO NOT EDIT -__version__ = "" +__version__ = "1000.0"