Fix some bugs in EIG calculator.

Now we can use either standard inversion or power method to
compute the spectral radii and plot the results.

src/msspec/calcio.py

@@ -594,7 +594,7 @@ class SpecIO(object):
         else:
             nlines = 1
         table = np.chararray((nlines, 4), unicode=True)
-        table[:] = str(p.get_parameter('eigval_ispectrum_ne').default)
+        table[:] = str(p.get_parameter('eigval_ispectrum_ne').value)
         for i in range(nlines):
             line = create_line("I_SPECTRUM(NE)")
             for j, o in enumerate((9, 19, 29, 39)):
@@ -944,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((2,-1))
         return data
 
     def load_facdif(self, filename='facdif1.dat'):

src/msspec/calculator.py

@@ -916,6 +916,10 @@ class _EIG(_MSCALCULATOR):
 
         self.spectroscopy_parameters.set_parameter('level', level)
 
+        # ensure that if the alogorithm is 'inversion' the IPWM must be set to 0
+        if self.global_parameters.algorithm.lower() == 'inversion':
+            self.spec_parameters.eigval_ipwm = 0
+
         self.get_tmatrix()
         self.run_spec()
 
@@ -939,7 +943,8 @@ class _EIG(_MSCALCULATOR):
         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 = dset.add_view("Spectral radius", title=title, xlabel=xlabel, ylabel=ylabel,
+                             marker='o')
         view.select('energy', 'spectral_radius')
 
         # add the cluster to the dataset

src/msspec/parameters.py

@@ -575,7 +575,7 @@ class SpecParameters(BaseParameters):
                       fmt='d'),
             Parameter('calc_igr', types=int, limits=[0, 2], default=0, fmt='d'),
 
-            Parameter('calc_iren', types=int, limits=[0, 4], default=1, fmt='d'),
+            Parameter('calc_iren', types=int, limits=[0, 5], default=1, fmt='d'),
             Parameter('calc_nren', types=int, limits=[1, None], default=1, fmt='d'),
             Parameter('calc_renr', types=float, limits=[None, None], default=1., fmt='.3f'),
             Parameter('calc_reni', types=float, limits=[None, None], default=0., fmt='.3f'),
@@ -1384,18 +1384,14 @@ class CalculationParameters(BaseParameters):
             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',
-                                                              'Z_n', 'K^2'), 
+                                                              'Z_n', 'Pi_1', 'Lowdin'), 
 		              types=(type(None), str), default=None,
                       doc=textwrap.dedent("""
             Enable the calculation of the coefficients for the renormalization of
             the multiple scattering series.
-            You can choose to renormalize in terms of the :math:`\\Sigma_n`, the 
+            You can choose to renormalize in terms of the :math:`G_n`, the 
-            :math:`G_n`, the :math:`Z_n` or the Lowdin :math:`K^2` matrices.
+            :math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` 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("""
@@ -1530,18 +1526,21 @@ class CalculationParameters(BaseParameters):
                     'G_n'    : 1,
                     'Sigma_n': 2,
                     'Z_n'    : 3,
-                    'K^2'    : 4}
+                    'Pi_1'   : 4,
+                    'Lowdin' : 5}
         # 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') )
+                    #assert( p.value in (None, 'Sigma_n', 'G_n') ) 
+                    assert( p.value in p.allowed_values)
             elif (self.global_parameters.spectroscopy == 'EIG'
-                  and self.global_parameters.algorithm == 'inversion'):
+                and self.global_parameters.algorithm == 'power'):
                     assert( p.value in p.allowed_values)
             else:
                 assert( p.value is None )
+            print(iren_map[p.value])
             self.spec_parameters.calc_iren = iren_map[p.value]
             LOGGER.info("Renormalization activated with \'{}\' method".format(p.value))
         except AssertionError:
@@ -1794,6 +1793,14 @@ class EIGParameters(BaseParameters):
             Parameter('kernel_matrix_spectrum', types=(bool,), default=False, doc=textwrap.dedent("""
             Whether to output the kernel matrix spectrum for each energy point.
             """)),
+            Parameter('method', types=(str,), default='EPSI',
+                      allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG',
+                                      'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER',
+                                      'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV',
+                                      'ELEV', 'EULE', 'GBWT', 'VARI', 'ITHE',
+                                      'EALG'],
+                      doc=textwrap.dedent("""The convergence acceleration scheme to be used.
+            """)),
         )
         BaseParameters.__init__(self)
         self.add_parameters(*parameters)
@@ -1811,3 +1818,12 @@ class EIGParameters(BaseParameters):
         self.spec_parameters.ped_li = li
         self.spec_parameters.ped_so = so
         self.spec_parameters.extra_level = p.value
+
+    def bind_method(self, p):
+        self.spec_parameters.eigval_method = p.value
+
+    def bind_kernel_matrix_spectrum(self, p):
+        value = int(p.value)
+        self.spec_parameters.eigval_ispectrum_ne = value
+
+        

src/msspec/spec/fortran/eig/pw/do_main.f

@@ -1229,7 +1229,8 @@ c          ENDIF
 C
       IF((ISOM.NE.0).OR.(NFICHLEC.EQ.1)) CLOSE(IUO1)
       IF(ISOM.NE.0) CLOSE(IUO2)
-      STOP
+CST      STOP
+      return
 C
    1  WRITE(IUO1,60)
       STOP