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

@@ -1,6 +1,10 @@
 __pycache__
 *.py[cod]
 *.so
+*.o
+*.mod
 *.bak
+*.target
+htmlcov
 package/
 src/msspec/results.txt

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)
+	@$(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
+#	remove *.lsm file and src_
-	@mv ./src/msspec/version.py.bak ./src/msspec/version.py
 	@rm ./msspec.lsm
+	@rm -rf ./src_
 
 
 version:

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.append(os.path.abspath('./'))
-sys.path.insert(0, os.path.abspath('../src'))
+sys.path.append(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('../'))
-
 
 import custom
 

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

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"
 

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:
+        #if len(data.shape) <= 1:
-            data = data.reshape((1, data.shape[0]))
+        #    data = data.reshape((1, data.shape[0]))
         return data
 
     def load_facdif(self, filename='facdif1.dat'):

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.10', 'cluster.clu')
-            os.rename('fort.35', 'tmatrix.tl')
+            #os.rename('fort.35', 'tmatrix.tl')
-            os.rename('fort.55', 'tmatrix.rad')
+            #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
 

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:

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
+            You can choose to renormalize in terms of the :math:`\\Sigma_n`, the 
-            renormalize in terms of the Sigma_n or G_n matrices.""")),
+            :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:
+        iren_map = {None     : 0,
-            self.spec_parameters.calc_iren = 0
+                    'G_n'    : 1,
-        else:
+                    'Sigma_n': 2,
-            if p.value.lower() == 'G_n'.lower():
+                    'Z_n'    : 3,
-                self.spec_parameters.calc_iren = 1
+                    'K^2'    : 4}
-            elif p.value.lower() == 'Sigma_n'.lower():
+        # Check that the method is neither 'Z_n' nor 'K^2' for other
-                self.spec_parameters.calc_iren = 2
+        # '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)

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

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
+    return cluster

src/msspec/version.py

@@ -1,3 +1,3 @@
 # This file is updated at package creation
 # DO NOT EDIT
-__version__ = ""
+__version__ = "1000.0"