From ff600a2e3db1fe9dccf3a4639e25753fd8dc4b9b Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Wed, 11 Mar 2020 11:08:17 +0100 Subject: [PATCH] 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 | 6 ++-- src/msspec/calculator.py | 7 ++++- src/msspec/parameters.py | 40 +++++++++++++++++------- src/msspec/spec/fortran/eig/pw/do_main.f | 3 +- 4 files changed, 39 insertions(+), 17 deletions(-) diff --git a/src/msspec/calcio.py b/src/msspec/calcio.py index a4c8004..ce14e09 100644 --- a/src/msspec/calcio.py +++ b/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: - # data = data.reshape((1, data.shape[0])) + if len(data.shape) <= 1: + data = data.reshape((2,-1)) return data def load_facdif(self, filename='facdif1.dat'): diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 4150cd1..33a0a1a 100644 --- a/src/msspec/calculator.py +++ b/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 diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 3b98f3b..96fa276 100644 --- a/src/msspec/parameters.py +++ b/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 - :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.""")), + You can choose to renormalize in terms of the :math:`G_n`, the + :math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin + :math:`K^2` matrices""")), 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 + + diff --git a/src/msspec/spec/fortran/eig/pw/do_main.f b/src/msspec/spec/fortran/eig/pw/do_main.f index 32d8d6f..c2c8bdf 100644 --- a/src/msspec/spec/fortran/eig/pw/do_main.f +++ b/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