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

@@ -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 _phd_se_noso_nosp_nosym, _phd_mi_noso_nosp_nosym
-from msspec.spec.fortran import eig_mi
+from msspec.spec.fortran import _eig_mi
-from msspec.spec.fortran import eig_pw
+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

src/msspec/misc.py

@@ -1,86 +1,88 @@
-# coding: utf-8
+# coding: utf-8
-"""
+# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : #
-Module misc
+"""
-===========
+Module misc
-
+===========
-"""
+
-import logging
+"""
-from pint import UnitRegistry
+import logging
-import numpy as np
+from pint import UnitRegistry
-import inspect
+import numpy as np
-import re
+import inspect
-
+import re
-class XRaySource(object):
+import os
-    MG_KALPHA = 1253.6
+
-    AL_KALPHA = 1486.6
+class XRaySource(object):
-    def __init__(self):
+    MG_KALPHA = 1253.6
-        pass
+    AL_KALPHA = 1486.6
-
+    def __init__(self):
-UREG = UnitRegistry()
+        pass
-UREG.define('rydberg = c * h * rydberg_constant = Ry')
+
-UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0')
+UREG = UnitRegistry()
-
+UREG.define('rydberg = c * h * rydberg_constant = Ry')
-logging.basicConfig(level=logging.INFO)
+UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0')
-LOGGER = logging.getLogger('msspec')
+
-
+logging.basicConfig(level=logging.INFO)
-np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5)
+LOGGER = logging.getLogger('msspec')
-
+
-def set_log_level(level):
+np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5)
-    lvl = getattr(logging, level.upper())
+
-    LOGGER.setLevel(lvl)
+def set_log_level(level):
-
+    lvl = getattr(logging, level.upper())
-def set_log_output(stream):
+    LOGGER.setLevel(lvl)
-    LOGGER.parent.handlers[0].stream = stream
+
-
+def set_log_output(stream):
-def get_call_info(frame):
+    LOGGER.parent.handlers[0].stream = stream
-    args, varargs, varkw, loc = inspect.getargvalues(frame)
+
-    _, _, function, _, _ = inspect.getframeinfo(frame)
+def get_call_info(frame):
-    s = '%s called with:\n' % function
+    args, varargs, varkw, loc = inspect.getargvalues(frame)
-
+    _, _, function, _, _ = inspect.getframeinfo(frame)
-    for kws in (args, varargs, varkw):
+    s = '%s called with:\n' % function
-        if kws != None:
+
-            if isinstance(kws, (tuple, list)):
+    for kws in (args, varargs, varkw):
-                for kw in kws:
+        if kws != None:
-                    s += '\t\t%s = %r\n' % (kw, loc[kw])
+            if isinstance(kws, (tuple, list)):
-            else:
+                for kw in kws:
-                s += '\t\t%s = %r\n' % (kws, loc[kws])
+                    s += '\t\t%s = %r\n' % (kw, loc[kw])
-    return s.strip()
+            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
+def log_process_output(process, logger=None, severity=logging.INFO):
-    else:
+    if logger == None:
-        logger = logging.getLogger(logger)
+        logger = logging
-        logger.setLevel(LOGGER.getEffectiveLevel())
+    else:
-
+        logger = logging.getLogger(logger)
-    for line in iter(process.stdout.readline, b''):
+        logger.setLevel(LOGGER.getEffectiveLevel())
-        logger.log(severity, line.rstrip('\n'))
+
-
+    for line in iter(process.stdout.readline, b''):
-    process.stdout.close()
+        logger.log(severity, line.rstrip('\n'))
-    process.wait()
+
-
+    process.stdout.close()
-
+    process.wait()
-def get_level_from_electron_configuration(notation):
+
-    l_letters = 'spdfghiklmnoqrtuv'
+
-    n_letters = 'klmnopqrstuv'
+def get_level_from_electron_configuration(notation):
-    pattern = re.compile(r'(?P<n>\d+)(?P<l>[%s%s])'
+    l_letters = 'spdfghiklmnoqrtuv'
-                         r'((?P<j>\d)/2)?' % (l_letters, l_letters.upper()))
+    n_letters = 'klmnopqrstuv'
-    m = pattern.match(notation)
+    pattern = re.compile(r'(?P<n>\d+)(?P<l>[%s%s])'
-    assert m, 'Not a valid notation'
+                         r'((?P<j>\d)/2)?' % (l_letters, l_letters.upper()))
-    n, l, _, j = m.groups()
+    m = pattern.match(notation)
-    n = int(n)
+    assert m, 'Not a valid notation'
-    l = l_letters.index(l.lower())
+    n, l, _, j = m.groups()
-    assert (l < n), 'Invalid l'
+    n = int(n)
-    j1 = abs(l - 0.5)
+    l = l_letters.index(l.lower())
-    j2 = abs(l + 0.5)
+    assert (l < n), 'Invalid l'
-    if j:
+    j1 = abs(l - 0.5)
-        j = int(j) / 2.
+    j2 = abs(l + 0.5)
-    else:
+    if j:
-        j = j1
+        j = int(j) / 2.
-    assert j in (j1, j2), 'Invalid j'
+    else:
-    letter = n_letters[n - 1]
+        j = j1
-    subscript = str(int((2 * (l + j) + 1) / 2))
+    assert j in (j1, j2), 'Invalid j'
-    if letter == 'k':
+    letter = n_letters[n - 1]
-        subscript = ''
+    subscript = str(int((2 * (l + j) + 1) / 2))
-    return '{}{}'.format(letter, subscript)
+    if letter == 'k':
+        subscript = ''
+    return '{}{}'.format(letter, subscript)

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'}
+        mapping = {'muffin_tin': 'in', 'lmto': 'ex', 'spkkr': 'ex', 'msf': 'ex'}
-        self.phagen_parameters.potgen = mapping[p.value]
+        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':
+        else:
-            self.phagen_parameters.potype = 'lmto'
+            self.phagen_parameters.potype = potential
+
+    def bind_potential_file(self, p):
+        pass 
 
     def bind_imaginery_part(self, p):
         self.phagen_parameters.gamma = p.value