From 927ac8a8a30a6b75e0ae50a5b47ad5dfa6063b23 Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Thu, 26 Mar 2020 18:26:27 +0100 Subject: [PATCH] Allow the use of external potential files. A new keyword in TMatrixParameters allows to use an external file for the potential energy of atoms. It should be used like this: calc.tmatrix_parameters.potential = 'spkkr' calc.tmatrix_parameters.potential_file = 'Cu.pot' # the name does not matter --- src/msspec/calculator.py | 33 +++++--- src/msspec/misc.py | 174 ++++++++++++++++++++------------------- src/msspec/parameters.py | 31 ++++--- 3 files changed, 130 insertions(+), 108 deletions(-) diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 33a0a1a..8d00198 100644 --- a/src/msspec/calculator.py +++ b/src/msspec/calculator.py @@ -1,5 +1,5 @@ # coding: utf-8 -# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a: +# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a tw=80: """ Module calculator ================= @@ -47,7 +47,8 @@ import numpy as np from msspec import iodata from msspec.data import electron_be from msspec.config import Config -from msspec.misc import (UREG, LOGGER, get_call_info, get_level_from_electron_configuration, +from msspec.misc import (UREG, LOGGER, get_call_info, + get_level_from_electron_configuration, XRaySource, set_log_output, log_process_output) from msspec.utils import get_atom_index @@ -62,9 +63,9 @@ from msspec.calcio import PhagenIO, SpecIO from msspec.phagen.fortran.libphagen import main as do_phagen -from msspec.spec.fortran import phd_se_noso_nosp_nosym, phd_mi_noso_nosp_nosym -from msspec.spec.fortran import eig_mi -from msspec.spec.fortran import eig_pw +from msspec.spec.fortran import _phd_se_noso_nosp_nosym, _phd_mi_noso_nosp_nosym +from msspec.spec.fortran import _eig_mi +from msspec.spec.fortran import _eig_pw from terminaltables.ascii_table import AsciiTable @@ -289,14 +290,21 @@ class _MSCALCULATOR(Calculator): # run phagen #self._make('tmatrix') os.chdir(os.path.join(self.tmp_folder, 'output')) + # copy external potential file to the workspace + if self.tmatrix_parameters.potential in ('lmto', 'spkkr', 'msf'): + fname = os.path.join(self.init_folder, + self.tmatrix_parameters.potential_file) + shutil.copy(fname, 'fort.2') 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') - shutil.copy('fort.10', 'cluster.clu') + #shutil.copy('fort.10', 'cluster.clu') + shutil.copy('clus/clus.out', 'cluster.clu') shutil.copy('fort.35', 'tmatrix.tl') shutil.copy('fort.55', 'tmatrix.rad') + os.chdir(self.init_folder) @@ -362,9 +370,9 @@ class _MSCALCULATOR(Calculator): # Get the spec function to launch if self.global_parameters.spectroscopy == 'PED': if self.global_parameters.algorithm == 'expansion': - do_spec = phd_se_noso_nosp_nosym.run + do_spec = _phd_se_noso_nosp_nosym.run elif self.global_parameters.algorithm == 'inversion': - do_spec = phd_mi_noso_nosp_nosym.run + do_spec = _phd_mi_noso_nosp_nosym.run else: LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not " "an allowed combination.".format(self.global_parameters.spectroscopy, @@ -372,9 +380,9 @@ class _MSCALCULATOR(Calculator): raise ValueError elif self.global_parameters.spectroscopy == 'EIG': if self.global_parameters.algorithm == 'inversion': - do_spec = eig_mi.run + do_spec = _eig_mi.run elif self.global_parameters.algorithm == 'power': - do_spec = eig_pw.run + do_spec = _eig_pw.run else: LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not " "an allowed combination.".format(self.global_parameters.spectroscopy, @@ -745,7 +753,10 @@ class _PED(_MSCALCULATOR): filename = os.path.join(self.tmp_folder, 'output/cluster.clu') self.phagenio.load_cluster_file(filename) - filename = os.path.join(self.tmp_folder, 'bin/plot/plot_vc.dat') + if self.phagen_parameters.potgen in ('in'): + filename = os.path.join(self.tmp_folder, 'output/plot/plot_vc.dat') + else: + filename = os.path.join(self.tmp_folder, 'output/plot/plot_v.dat') pot_data = self.phagenio.load_potential_file(filename) cluster = self.phagen_parameters.get_parameter('atoms').value diff --git a/src/msspec/misc.py b/src/msspec/misc.py index b828b47..6054210 100644 --- a/src/msspec/misc.py +++ b/src/msspec/misc.py @@ -1,86 +1,88 @@ -# coding: utf-8 -""" -Module misc -=========== - -""" -import logging -from pint import UnitRegistry -import numpy as np -import inspect -import re - -class XRaySource(object): - MG_KALPHA = 1253.6 - AL_KALPHA = 1486.6 - def __init__(self): - pass - -UREG = UnitRegistry() -UREG.define('rydberg = c * h * rydberg_constant = Ry') -UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0') - -logging.basicConfig(level=logging.INFO) -LOGGER = logging.getLogger('msspec') - -np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5) - -def set_log_level(level): - lvl = getattr(logging, level.upper()) - LOGGER.setLevel(lvl) - -def set_log_output(stream): - LOGGER.parent.handlers[0].stream = stream - -def get_call_info(frame): - args, varargs, varkw, loc = inspect.getargvalues(frame) - _, _, function, _, _ = inspect.getframeinfo(frame) - s = '%s called with:\n' % function - - for kws in (args, varargs, varkw): - if kws != None: - if isinstance(kws, (tuple, list)): - for kw in kws: - s += '\t\t%s = %r\n' % (kw, loc[kw]) - else: - s += '\t\t%s = %r\n' % (kws, loc[kws]) - return s.strip() - - -def log_process_output(process, logger=None, severity=logging.INFO): - if logger == None: - logger = logging - else: - logger = logging.getLogger(logger) - logger.setLevel(LOGGER.getEffectiveLevel()) - - for line in iter(process.stdout.readline, b''): - logger.log(severity, line.rstrip('\n')) - - process.stdout.close() - process.wait() - - -def get_level_from_electron_configuration(notation): - l_letters = 'spdfghiklmnoqrtuv' - n_letters = 'klmnopqrstuv' - pattern = re.compile(r'(?P\d+)(?P[%s%s])' - r'((?P\d)/2)?' % (l_letters, l_letters.upper())) - m = pattern.match(notation) - assert m, 'Not a valid notation' - n, l, _, j = m.groups() - n = int(n) - l = l_letters.index(l.lower()) - assert (l < n), 'Invalid l' - j1 = abs(l - 0.5) - j2 = abs(l + 0.5) - if j: - j = int(j) / 2. - else: - j = j1 - assert j in (j1, j2), 'Invalid j' - letter = n_letters[n - 1] - subscript = str(int((2 * (l + j) + 1) / 2)) - if letter == 'k': - subscript = '' - return '{}{}'.format(letter, subscript) +# coding: utf-8 +# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : # +""" +Module misc +=========== + +""" +import logging +from pint import UnitRegistry +import numpy as np +import inspect +import re +import os + +class XRaySource(object): + MG_KALPHA = 1253.6 + AL_KALPHA = 1486.6 + def __init__(self): + pass + +UREG = UnitRegistry() +UREG.define('rydberg = c * h * rydberg_constant = Ry') +UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0') + +logging.basicConfig(level=logging.INFO) +LOGGER = logging.getLogger('msspec') + +np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5) + +def set_log_level(level): + lvl = getattr(logging, level.upper()) + LOGGER.setLevel(lvl) + +def set_log_output(stream): + LOGGER.parent.handlers[0].stream = stream + +def get_call_info(frame): + args, varargs, varkw, loc = inspect.getargvalues(frame) + _, _, function, _, _ = inspect.getframeinfo(frame) + s = '%s called with:\n' % function + + for kws in (args, varargs, varkw): + if kws != None: + if isinstance(kws, (tuple, list)): + for kw in kws: + s += '\t\t%s = %r\n' % (kw, loc[kw]) + else: + s += '\t\t%s = %r\n' % (kws, loc[kws]) + return s.strip() + + +def log_process_output(process, logger=None, severity=logging.INFO): + if logger == None: + logger = logging + else: + logger = logging.getLogger(logger) + logger.setLevel(LOGGER.getEffectiveLevel()) + + for line in iter(process.stdout.readline, b''): + logger.log(severity, line.rstrip('\n')) + + process.stdout.close() + process.wait() + + +def get_level_from_electron_configuration(notation): + l_letters = 'spdfghiklmnoqrtuv' + n_letters = 'klmnopqrstuv' + pattern = re.compile(r'(?P\d+)(?P[%s%s])' + r'((?P\d)/2)?' % (l_letters, l_letters.upper())) + m = pattern.match(notation) + assert m, 'Not a valid notation' + n, l, _, j = m.groups() + n = int(n) + l = l_letters.index(l.lower()) + assert (l < n), 'Invalid l' + j1 = abs(l - 0.5) + j2 = abs(l + 0.5) + if j: + j = int(j) / 2. + else: + j = j1 + assert j in (j1, j2), 'Invalid j' + letter = n_letters[n - 1] + subscript = str(int((2 * (l + j) + 1) / 2)) + if letter == 'k': + subscript = '' + return '{}{}'.format(letter, subscript) diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index d848c56..fea0964 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -265,6 +265,7 @@ class BaseParameters(object): class PhagenParameters(BaseParameters): def __init__(self): parameters = ( + Parameter('version', types=(str,), fmt='>3s', default='1.1'), Parameter('calctype', allowed_values=('xpd', 'xas', 'aed', 'led', 'rex', 'els', 'e2e'#, 'dos' ), @@ -273,7 +274,7 @@ class PhagenParameters(BaseParameters): types=(str,), default='cis'), Parameter('coor', allowed_values=('angs', 'au'), types=(str,), default='angs'), - Parameter('enunit', allowed_values=('ryd', 'eV'), fmt='3>s', + Parameter('enunit', allowed_values=('ryd', 'eV'), fmt='>3s', types=(str,), default='ryd'), Parameter('einc', types=(int, float), limits=(0, None), fmt='.1f', default=700.), @@ -281,7 +282,7 @@ class PhagenParameters(BaseParameters): default=580.), Parameter('scangl', types=(int, float), limits=(0, 360), fmt='.2f', default=0.), - Parameter('lambda', types=(int, float), fmt='.2f', default=20.), + Parameter('lambda', types=(int, float), fmt='.1f', default=0.), Parameter('emin', types=(int, float), limits=(0, None), fmt='.4f', default=13.5236), Parameter('emax', types=(int, float), limits=(0, None), fmt='.4f', @@ -294,9 +295,9 @@ class PhagenParameters(BaseParameters): default='in'), Parameter('potype', allowed_values=('hdrel', 'hedin', 'xalph', 'dhrel', 'dhcmp', - #'lmto', 'msf', 'spkkr' + 'lmto', 'msf', 'spkkr' ), - types=(str,), fmt='5>s', default='hedin'), + types=(str,), fmt='>5s', default='hedin'), Parameter('relc', allowed_values=('nr', 'sr', 'so'), types=(str,), default='nr'), Parameter('norman', allowed_values=('stdcrm', 'scaled', 'extrad'), @@ -308,9 +309,9 @@ class PhagenParameters(BaseParameters): Parameter('charelx', allowed_values=('ex', 'gs'), types=(str,), default='gs'), Parameter('ionzst', allowed_values=('neutral', 'ionic'), - types=(str,), fmt='7>s', default='neutral'), + types=(str,), fmt='>7s', default='neutral'), Parameter('eikappr', allowed_values=('yes', 'no'), types=(str,), - fmt='3>s', default='no'), + fmt='>3s', default='no'), Parameter('db', types=(int, float), fmt='.2f', default=0.01), Parameter('optrsh', allowed_values=('y', 'n'), types=(str,), default='n'), @@ -885,11 +886,12 @@ class TMatrixParameters(BaseParameters): def __init__(self, phagen_parameters): parameters = ( Parameter('potential', types=str, - allowed_values=('muffin_tin', 'lmto'), + allowed_values=('muffin_tin', 'lmto', 'spkkr', 'msf'), default='muffin_tin', doc=textwrap.dedent(""" This option allows to choose which kind of potential is used to compute the T-Matrix. For now, only the internal *Muffin-Tin* potential is supported. """)), + Parameter('potential_file', types=(str,type(None)), default=None), Parameter('exchange_correlation', types=str, allowed_values=('hedin_lundqvist_real', 'hedin_lundqvist_complex', @@ -955,8 +957,12 @@ class TMatrixParameters(BaseParameters): self.freeze() def bind_potential(self, p): - mapping = {'muffin_tin': 'in', 'lmto': 'ex'} - self.phagen_parameters.potgen = mapping[p.value] + mapping = {'muffin_tin': 'in', 'lmto': 'ex', 'spkkr': 'ex', 'msf': 'ex'} + if p.value.lower() in ('muffin_tin'): + self.phagen_parameters.potgen = 'in' + elif p.value.lower() in ('spkkr', 'lmto', 'msf'): + self.phagen_parameters.potgen = 'ex' + self.phagen_parameters.potype = p.value.lower() def bind_exchange_correlation(self, p): potential = self.get_parameter('potential').value @@ -969,8 +975,11 @@ class TMatrixParameters(BaseParameters): 'dirac_hara_complex': 'dhcmp' } self.phagen_parameters.potype = mapping[p.value] - elif potential == 'lmto': - self.phagen_parameters.potype = 'lmto' + else: + self.phagen_parameters.potype = potential + + def bind_potential_file(self, p): + pass def bind_imaginery_part(self, p): self.phagen_parameters.gamma = p.value