epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Add support for SPRKKR potential. 

This is a first version for this option. Some work has still to be
done...
This commit is contained in:
Sylvain Tricot 2020-04-10 17:36:25 +02:00
parent 81459cfcfc
commit 8ebfd624e1
13 changed files with 1278 additions and 1259 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

302
doc/source/tutorials/copper/log.log
View File

@ -1,302 +0,0 @@
INFO:msspec:Getting the TMatrix...
_________________________________________________________________
PHAGEN
_________________________________________________________________
parameters for this xpd calculation:
-----------------------------------------------------------------
edge= k
potype= hedin norman= stdcrm absorber= 1
coor= angs emin= 13.52 Ry emax= 13.52 Ry
delta= 0.300 Ry gamma= 0.03 Ry eftri= 0.000 Ry
ionization state : neutral
final state potential generated internally
Computes the T-matrix and radial matrix elements
coordinates in angstroms Radii
-----------------------------------------------------------------
Cu 29 0.0000 0.0000 -3.6000 0.0000 0.0000
Cu 29 -1.8000 -1.8000 0.0000 0.0000 0.0000
Cu 29 -1.8000 0.0000 -1.8000 0.0000 0.0000
Cu 29 -3.6000 0.0000 0.0000 0.0000 0.0000
Cu 29 -1.8000 1.8000 0.0000 0.0000 0.0000
Cu 29 0.0000 -1.8000 -1.8000 0.0000 0.0000
Cu 29 0.0000 -3.6000 0.0000 0.0000 0.0000
Cu 29 1.8000 -1.8000 0.0000 0.0000 0.0000
Cu 29 0.0000 1.8000 -1.8000 0.0000 0.0000
Cu 29 1.8000 0.0000 -1.8000 0.0000 0.0000
Cu 29 0.0000 0.0000 0.0000 0.0000 0.0000
Cu 29 1.8000 1.8000 0.0000 0.0000 0.0000
Cu 29 0.0000 3.6000 0.0000 0.0000 0.0000
Cu 29 3.6000 0.0000 0.0000 0.0000 0.0000
-----------------------------------------------------------------
** enter calphas **
---
total energy for atom in ground state
total energy for atom with a hole in k edge
calculated ionization energy for edge k = 8994.86914 eV
---
calculated ionization potential (ryd) = 661.387451
---
symmetrizing coordinates...
symmetrized atomic coordinates of cluster
position
atom no. x y z eq
1 osph 0 0.0000 0.0000 0.0000 0
2 Cu 29 0.0000 0.0000 -5.3452 0
3 Cu 29 -3.4015 -3.4015 1.4578 0
4 Cu 29 -3.4015 0.0000 -1.9437 0
5 Cu 29 -6.8030 0.0000 1.4578 0
6 Cu 29 0.0000 0.0000 1.4578 0
7 Cu 29 -3.4015 3.4015 1.4578 3
8 Cu 29 3.4015 -3.4015 1.4578 3
9 Cu 29 3.4015 3.4015 1.4578 3
10 Cu 29 0.0000 -3.4015 -1.9437 4
11 Cu 29 0.0000 3.4015 -1.9437 4
12 Cu 29 3.4015 0.0000 -1.9437 4
13 Cu 29 0.0000 -6.8030 1.4578 5
14 Cu 29 0.0000 6.8030 1.4578 5
15 Cu 29 6.8030 0.0000 1.4578 5
computing muffin tin potential and phase shifts
generating core state wavefunction
generating final potential (complex hedin-lundqvist exchange)
MT radii for Hydrogen atoms determined by stdcrm unless other options are specified
-----------------------------------------------------------------
i rs(i) i=1,natoms
1 9.28 2 2.50 3 2.46 4 2.35 5 2.48 6 2.33 7 2.46 8 2.46
9 2.46 10 2.35 11 2.35 12 2.35 13 2.48 14 2.48 15 2.48
N.B.: Order of atoms as reshuffled by symmetry routines
-----------------------------------------------------------------
input value for coulomb interst. potential = -0.699999988
and interstitial rs = 3.00000000
lower bound for coulomb interst. potential = -0.290025890
and for interst. rs = 2.54297900
lmax assignment based on l_max = r_mt * k_e + 2
where e is the running energy
optimal lmax chosen according to the running energy e for each atom
number of centers= 14
starting potentials and/or charge densities written to file 13
symmetry information generated internally
symmetry information written to file 14
core initial state of type: 1s1/2
fermi level = -0.17690
generating t_l (for030) and atomic cross section (for050)
gamma = 0.030000 rsint = 3.944844
check in subroutine cont
order of neighb. -- symb. -- dist. from absorber
3 Cu 4.81045771
5 Cu 6.80301476
2 Cu 8.33195782
4 Cu 9.62091541
-----------------------------------------------------------------
1 Cu 0.000000
3 Cu 4.810458
5 Cu 6.803015
2 Cu 8.331958
4 Cu 9.620915
1 Cu 0.000000
3 Cu 4.810458
5 Cu 6.803015
2 Cu 8.331958
4 Cu 9.620915
irho = 2 entering vxc to calculate energy dependent exchange
energy dependent vcon = (-0.181322575,0.183172509) at energy 13.5235996
calculating atomic t-matrix elements atm(n)
check orthogonality between core and continuum state
scalar product between core and continuum state = (0.215178907,-7.336383685E-03)
sqrt(xe/pi) = (1.08555424,-3.627028549E-03)
check ionization potential: 661.387451
value of the mean free path:
-----------------------------------------------------------------
average mean free path due to finite gamma: mfp = 8.81667 angstrom at energy 13.52360
average mean free path in the cluster : mfp = 8.71420 angstrom at energy 13.52360
total mean free path due to Im V and gamma: mfp = 4.38257 angstrom at energy 13.52360
-----------------------------------------------------------------
** phagen terminated normally **
INFO:msspec:Getting the TMatrix...
_________________________________________________________________
PHAGEN
_________________________________________________________________
parameters for this xpd calculation:
-----------------------------------------------------------------
edge= k
potype= hedin norman= stdcrm absorber= 1
coor= angs emin= 13.52 Ry emax= 13.52 Ry
delta= 0.300 Ry gamma= 0.03 Ry eftri= 0.000 Ry
ionization state : neutral
final state potential generated internally
Computes the T-matrix and radial matrix elements
coordinates in angstroms Radii
-----------------------------------------------------------------
Cu 29 0.0000 0.0000 -3.8000 0.0000 0.0000
Cu 29 -1.9000 -1.9000 0.0000 0.0000 0.0000
Cu 29 -1.9000 0.0000 -1.9000 0.0000 0.0000
Cu 29 -3.8000 0.0000 0.0000 0.0000 0.0000
Cu 29 -1.9000 1.9000 0.0000 0.0000 0.0000
Cu 29 0.0000 -1.9000 -1.9000 0.0000 0.0000
Cu 29 0.0000 -3.8000 0.0000 0.0000 0.0000
Cu 29 1.9000 -1.9000 0.0000 0.0000 0.0000
Cu 29 0.0000 1.9000 -1.9000 0.0000 0.0000
Cu 29 1.9000 0.0000 -1.9000 0.0000 0.0000
Cu 29 0.0000 0.0000 0.0000 0.0000 0.0000
Cu 29 1.9000 1.9000 0.0000 0.0000 0.0000
Cu 29 0.0000 3.8000 0.0000 0.0000 0.0000
Cu 29 3.8000 0.0000 0.0000 0.0000 0.0000
-----------------------------------------------------------------
** enter calphas **
---
total energy for atom in ground state
total energy for atom with a hole in k edge
calculated ionization energy for edge k = 8994.86914 eV
---
calculated ionization potential (ryd) = 661.387451
---
symmetrizing coordinates...
symmetrized atomic coordinates of cluster
position
atom no. x y z eq
1 osph 0 0.0000 0.0000 0.0000 0
2 Cu 29 0.0000 0.0000 -5.6422 0
3 Cu 29 -3.5905 -3.5905 1.5388 0
4 Cu 29 -3.5905 0.0000 -2.0517 0
5 Cu 29 -7.1810 0.0000 1.5388 0
6 Cu 29 0.0000 0.0000 1.5388 0
7 Cu 29 -3.5905 3.5905 1.5388 3
8 Cu 29 3.5905 -3.5905 1.5388 3
9 Cu 29 3.5905 3.5905 1.5388 3
10 Cu 29 0.0000 -3.5905 -2.0517 4
11 Cu 29 0.0000 3.5905 -2.0517 4
12 Cu 29 3.5905 0.0000 -2.0517 4
13 Cu 29 0.0000 -7.1810 1.5388 5
14 Cu 29 0.0000 7.1810 1.5388 5
15 Cu 29 7.1810 0.0000 1.5388 5
computing muffin tin potential and phase shifts
generating core state wavefunction
generating final potential (complex hedin-lundqvist exchange)
MT radii for Hydrogen atoms determined by stdcrm unless other options are specified
-----------------------------------------------------------------
i rs(i) i=1,natoms
1 9.80 2 2.64 3 2.60 4 2.44 5 2.61 6 2.46 7 2.60 8 2.60
9 2.60 10 2.44 11 2.44 12 2.44 13 2.61 14 2.61 15 2.61
N.B.: Order of atoms as reshuffled by symmetry routines
-----------------------------------------------------------------
input value for coulomb interst. potential = -0.699999988
and interstitial rs = 3.00000000
lower bound for coulomb interst. potential = -0.232031077
and for interst. rs = 2.76609278
lmax assignment based on l_max = r_mt * k_e + 2
where e is the running energy
optimal lmax chosen according to the running energy e for each atom
number of centers= 14
starting potentials and/or charge densities written to file 13
symmetry information generated internally
symmetry information written to file 14
core initial state of type: 1s1/2
fermi level = -0.16924
generating t_l (for030) and atomic cross section (for050)
gamma = 0.030000 rsint = 4.292365
check in subroutine cont
order of neighb. -- symb. -- dist. from absorber
3 Cu 5.07770586
5 Cu 7.18096018
2 Cu 8.79484367
4 Cu 10.1554117
-----------------------------------------------------------------
1 Cu 0.000000
3 Cu 5.077706
5 Cu 7.180960
2 Cu 8.794844
4 Cu 10.155412
1 Cu 0.000000
3 Cu 5.077706
5 Cu 7.180960
2 Cu 8.794844
4 Cu 10.155412
irho = 2 entering vxc to calculate energy dependent exchange
energy dependent vcon = (-0.151319593,0.168841615) at energy 13.5235996
calculating atomic t-matrix elements atm(n)
check orthogonality between core and continuum state
scalar product between core and continuum state = (0.217915282,-9.193832986E-03)
sqrt(xe/pi) = (1.08495688,-3.348779166E-03)
check ionization potential: 661.387451
value of the mean free path:
-----------------------------------------------------------------
average mean free path due to finite gamma: mfp = 8.81667 angstrom at energy 13.52360
average mean free path in the cluster : mfp = 9.13179 angstrom at energy 13.52360
total mean free path due to Im V and gamma: mfp = 4.48573 angstrom at energy 13.52360
-----------------------------------------------------------------
** phagen terminated normally **

30
src/msspec/__init__.py
View File

@ -1,8 +1,30 @@
import ase #!/usr/bin/env python
from .version import __version__ #
from .misc import LOGGER # Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/__init__.py
# Last modified: ven. 10 avril 2020 17:22:12
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
import ase
from msspec.misc import LOGGER
from msspec.version import __version__
__sha__ = "$Id$"
def init_msspec(): def init_msspec():
LOGGER.debug('Initialization of the msspec module') LOGGER.debug('Initialization of the msspec module')

145
src/msspec/calcio.py
View File

@ -1,4 +1,5 @@
# coding: utf-8 # coding: utf-8
# vim: set et ts=4 sw=4 fdm=indent mouse=a cc=+1 tw=80:
""" """
Module calcio Module calcio
@ -6,13 +7,16 @@ Module calcio
""" """
import numpy as np
from datetime import datetime from datetime import datetime
import ase
import re import re
import numpy as np
import ase
from ase import Atom, Atoms
from ase.data import chemical_symbols
from msspec.misc import UREG, LOGGER from msspec.misc import UREG, LOGGER
class PhagenIO(object): class PhagenIO(object):
def __init__(self, phagen_parameters, malloc_parameters): def __init__(self, phagen_parameters, malloc_parameters):
self.parameters = phagen_parameters self.parameters = phagen_parameters
@ -27,8 +31,9 @@ class PhagenIO(object):
self.nlmax = None self.nlmax = None
self.energies = None self.energies = None
def write_input_file(self, filename='input.ms'): def write_input_file(self, filename='input.ms'): # noqa: C901
# in input folder # in input folder
LOGGER.debug("Writing Phagen input file")
atoms = self.parameters.atoms atoms = self.parameters.atoms
if not atoms.has('mt_radii'): if not atoms.has('mt_radii'):
for a in atoms: for a in atoms:
@ -58,7 +63,7 @@ class PhagenIO(object):
# other sections # other sections
continue continue
value = parameter.value value = parameter.value
if isinstance(value, str) and not re.match('\..*\.', value): if isinstance(value, str) and not re.match(r'\..*\.', value):
s = ' {}=\'{:s}\''.format(name, str(parameter)) s = ' {}=\'{:s}\''.format(name, str(parameter))
else: else:
s = ' {}={:s}'.format(name, str(parameter)) s = ' {}={:s}'.format(name, str(parameter))
@ -167,7 +172,8 @@ class PhagenIO(object):
lines = pattern.split(content) lines = pattern.split(content)
# get the number of atoms (nat) and the number of energy points (ne) # get the number of atoms (nat) and the number of energy points (ne)
nat, ne, _, ipot, lmax_mode = list(map(int, content.split('\n')[0].split())) nat, ne, _, ipot, lmax_mode = list(map(int,
content.split('\n')[0].split()))
self.nat = nat self.nat = nat
self.ne = ne self.ne = ne
self.ipot = ipot self.ipot = ipot
@ -181,7 +187,7 @@ class PhagenIO(object):
if not re.match(r'^\d+$', n): if not re.match(r'^\d+$', n):
numbers.append(float(n)) numbers.append(float(n))
if len(numbers) > 0: if len(numbers) > 0:
array = np.array(numbers).reshape((-1, 4)) # pylint: disable=no-member array = np.array(numbers).reshape((-1, 4))
atom_data.append(array) atom_data.append(array)
# construct the data array # construct the data array
@ -265,8 +271,27 @@ class PhagenIO(object):
proto_index = int(data[i, -1]) proto_index = int(data[i, -1])
proto_indices.append(proto_index) proto_indices.append(proto_index)
atoms.set_array('proto_indices', np.array(proto_indices)) atoms.set_array('proto_indices', np.array(proto_indices))
self.nateqm = int(np.max([len(np.where(data[:,-1]==i)[0]) for i in range( self.nateqm = int(np.max([
1, self.nat + 1)])) len(np.where(data[:, -1] == i)[0]) for i in range(1, self.nat + 1)
]))
def load_prototypical_atoms(self, filename='molinpot3.out'):
with open(filename, 'r') as fd:
content = fd.readlines()
content = content[-len(self.parameters.atoms):]
cluster = Atoms()
for line in content:
datatxt = np.array(re.split(r'\s+', line.strip(' \n')))
neq = int(datatxt[-1])
if neq == 0:
Z = int(datatxt[1])
xyz = (datatxt[[2, 3, 4]].astype(float) * UREG.a0).to(
'angstrom').magnitude
sym = chemical_symbols[Z]
cluster += Atom(sym, position=xyz)
self.prototypical_atoms = cluster
return cluster
def load_potential_file(self, filename='plot_vc.dat'): def load_potential_file(self, filename='plot_vc.dat'):
a_index = 0 a_index = 0
@ -279,11 +304,13 @@ class PhagenIO(object):
a_index += 1 a_index += 1
d = d.split() d = d.split()
a = {'Symbol': d[1], 'distance': float(d[4]), a = {'Symbol': d[1], 'distance': float(d[4]),
'coord': np.array([float(d[7]), float(d[8]), float(d[9])]), 'coord': np.array([float(d[7]), float(d[8]),
float(d[9])]),
'index': int(a_index), 'values': []} 'index': int(a_index), 'values': []}
pot_data.append(a) pot_data.append(a)
else: else:
pot_data[a_index - 1]['values'].append(tuple(float(_) for _ in d.split())) pot_data[a_index - 1]['values'].append(
tuple(float(_) for _ in d.split()))
# convert the values list to a numpy array # convert the values list to a numpy array
for _pot_data in pot_data: for _pot_data in pot_data:
@ -292,27 +319,30 @@ class PhagenIO(object):
return pot_data return pot_data
class SpecIO(object): class SpecIO(object):
def __init__(self, parameters, malloc_parameters, phagenio): def __init__(self, parameters, malloc_parameters, phagenio):
self.parameters = parameters self.parameters = parameters
self.malloc_parameters = malloc_parameters self.malloc_parameters = malloc_parameters
self.phagenio = phagenio self.phagenio = phagenio
def write_input_file(self, filename='spec.dat'): def write_input_file(self, filename='spec.dat'): # noqa: C901
def title(t, shift=4, width=79, center=True): def title(t, shift=4, width=79, center=True):
if center: if center:
s = ('{}*{:^%ds}*\n' % (width - shift - 2)).format(' ' * shift, t) s = ('{}*{:^%ds}*\n' % (width - shift - 2)
).format(' ' * shift, t)
else: else:
s = ('{}*{:%ds}*\n' % (width - shift - 2)).format(' ' * shift, t) s = ('{}*{:%ds}*\n' % (width - shift - 2)
).format(' ' * shift, t)
return s return s
def rule(tabs=(5, 10, 10, 10, 10), symbol='=', shift=4, width=79): def rule(tabs=(5, 10, 10, 10, 10), symbol='=', shift=4, width=79):
s = ' ' * shift + '*' + symbol * (width - shift - 2) + '*\n' s = ' ' * shift + '*' + symbol * (width - shift - 2) + '*\n'
t = np.cumsum(tabs) + shift t = np.cumsum(tabs) + shift
l = list(s) sep = list(s)
for i in t: for i in t:
l[i] = '+' sep[i] = '+'
return ''.join(l) return ''.join(sep)
def fillstr(a, b, index, justify='left'): def fillstr(a, b, index, justify='left'):
alist = list(a) alist = list(a)
@ -426,7 +456,8 @@ class SpecIO(object):
content += rule() content += rule()
line = create_line("SPECTRO,ISPIN,IDICHR,IPOL") line = create_line("SPECTRO,ISPIN,IDICHR,IPOL")
line = fillstr(line, str(p.calctype_spectro), 9, 'left') line = fillstr(line, str(p.calctype_spectro), 9, 'left')
line = fillstr(line, str(p.get_parameter('calctype_ispin')), 19, 'left') line = fillstr(line,
str(p.get_parameter('calctype_ispin')), 19, 'left')
line = fillstr(line, str(p.calctype_idichr), 29, 'left') line = fillstr(line, str(p.calctype_idichr), 29, 'left')
line = fillstr(line, str(p.calctype_ipol), 39, 'left') line = fillstr(line, str(p.calctype_ipol), 39, 'left')
content += line content += line
@ -476,11 +507,13 @@ class SpecIO(object):
line = fillstr(line, str(p.get_parameter('ped_imod')), 9, 'decimal') line = fillstr(line, str(p.get_parameter('ped_imod')), 9, 'decimal')
line = fillstr(line, str(p.get_parameter('ped_imoy')), 19, 'decimal') line = fillstr(line, str(p.get_parameter('ped_imoy')), 19, 'decimal')
line = fillstr(line, str(p.get_parameter('ped_accept')), 29, 'decimal') line = fillstr(line, str(p.get_parameter('ped_accept')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('ped_ichkdir')), 39, 'decimal') line = fillstr(line,
str(p.get_parameter('ped_ichkdir')), 39, 'decimal')
content += line content += line
content += rule() content += rule()
content += title(' ' * 22 + 'LEED EXPERIMENTAL PARAMETERS :', center=False) content += title(' ' * 22 + 'LEED EXPERIMENTAL PARAMETERS :',
center=False)
content += rule() content += rule()
line = create_line("IPHI,ITHETA,IE,IFTHET") line = create_line("IPHI,ITHETA,IE,IFTHET")
line = fillstr(line, str(p.get_parameter('leed_iphi')), 9, 'left') line = fillstr(line, str(p.get_parameter('leed_iphi')), 9, 'left')
@ -496,48 +529,57 @@ class SpecIO(object):
content += line content += line
line = create_line("PHI0,THETA0,E0,R0") line = create_line("PHI0,THETA0,E0,R0")
line = fillstr(line, str(p.get_parameter('leed_phi0')), 9, 'decimal') line = fillstr(line, str(p.get_parameter('leed_phi0')), 9, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_theta0')), 19, 'decimal') line = fillstr(line,
str(p.get_parameter('leed_theta0')), 19, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_e0')), 29, 'decimal') line = fillstr(line, str(p.get_parameter('leed_e0')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_r0')), 39, 'decimal') line = fillstr(line, str(p.get_parameter('leed_r0')), 39, 'decimal')
content += line content += line
line = create_line("PHI1,THETA1,E1,R1") line = create_line("PHI1,THETA1,E1,R1")
line = fillstr(line, str(p.get_parameter('leed_phi1')), 9, 'decimal') line = fillstr(line, str(p.get_parameter('leed_phi1')), 9, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_theta1')), 19, 'decimal') line = fillstr(line,
str(p.get_parameter('leed_theta1')), 19, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_e1')), 29, 'decimal') line = fillstr(line, str(p.get_parameter('leed_e1')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_r1')), 39, 'decimal') line = fillstr(line, str(p.get_parameter('leed_r1')), 39, 'decimal')
content += line content += line
line = create_line("TH_INI,PHI_INI") line = create_line("TH_INI,PHI_INI")
line = fillstr(line, str(p.get_parameter('leed_thini')), 9, 'decimal') line = fillstr(line, str(p.get_parameter('leed_thini')), 9, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_phiini')), 19, 'decimal') line = fillstr(line,
str(p.get_parameter('leed_phiini')), 19, 'decimal')
content += line content += line
line = create_line("IMOD,IMOY,ACCEPT,ICHKDIR") line = create_line("IMOD,IMOY,ACCEPT,ICHKDIR")
line = fillstr(line, str(p.get_parameter('leed_imod')), 9, 'decimal') line = fillstr(line, str(p.get_parameter('leed_imod')), 9, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_imoy')), 19, 'decimal') line = fillstr(line, str(p.get_parameter('leed_imoy')), 19, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_accept')), 29, 'decimal') line = fillstr(line,
str(p.get_parameter('leed_accept')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('leed_ichkdir')), 39, line = fillstr(line, str(p.get_parameter('leed_ichkdir')), 39,
'decimal') 'decimal')
content += line content += line
content += rule() content += rule()
content += title(' ' * 21 + 'EXAFS EXPERIMENTAL PARAMETERS :', center=False) content += title(' ' * 21 + 'EXAFS EXPERIMENTAL PARAMETERS :',
center=False)
content += rule() content += rule()
line = create_line("EDGE,INITL,THLUM,PHILUM") line = create_line("EDGE,INITL,THLUM,PHILUM")
line = fillstr(line, str(p.get_parameter('exafs_edge')), 9, 'left') line = fillstr(line, str(p.get_parameter('exafs_edge')), 9, 'left')
line = fillstr(line, str(p.get_parameter('exafs_initl')), 19, 'left') line = fillstr(line, str(p.get_parameter('exafs_initl')), 19, 'left')
line = fillstr(line, str(p.get_parameter('exafs_thlum')), 29, 'decimal') line = fillstr(line,
str(p.get_parameter('exafs_thlum')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('exafs_philum')), 39, line = fillstr(line, str(p.get_parameter('exafs_philum')), 39,
'decimal') 'decimal')
content += line content += line
line = create_line("NE,EK_INI,EK_FIN,EPH_INI") line = create_line("NE,EK_INI,EK_FIN,EPH_INI")
line = fillstr(line, str(p.get_parameter('exafs_ne')), 9, 'left') line = fillstr(line, str(p.get_parameter('exafs_ne')), 9, 'left')
line = fillstr(line, str(p.get_parameter('exafs_ekini')), 19, 'decimal') line = fillstr(line,
line = fillstr(line, str(p.get_parameter('exafs_ekfin')), 29, 'decimal') str(p.get_parameter('exafs_ekini')), 19, 'decimal')
line = fillstr(line,
str(p.get_parameter('exafs_ekfin')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('exafs_ephini')), 39, line = fillstr(line, str(p.get_parameter('exafs_ephini')), 39,
'decimal') 'decimal')
content += line content += line
content += rule() content += rule()
content += title(' ' * 22 + 'AED EXPERIMENTAL PARAMETERS :', center=False) content += title(' ' * 22 + 'AED EXPERIMENTAL PARAMETERS :',
center=False)
content += rule() content += rule()
line = create_line("EDGE_C,EDGE_I,EDGE_A") line = create_line("EDGE_C,EDGE_I,EDGE_A")
line = fillstr(line, str(p.get_parameter('aed_edgec')), 9, 'left') line = fillstr(line, str(p.get_parameter('aed_edgec')), 9, 'left')
@ -605,7 +647,8 @@ class SpecIO(object):
line = fillstr(line, str(p.get_parameter('eigval_ipwm')), 9, 'left') line = fillstr(line, str(p.get_parameter('eigval_ipwm')), 9, 'left')
line = fillstr(line, str(p.get_parameter('eigval_method')), 19, 'left') line = fillstr(line, str(p.get_parameter('eigval_method')), 19, 'left')
line = fillstr(line, str(p.get_parameter('eigval_acc')), 29, 'decimal') line = fillstr(line, str(p.get_parameter('eigval_acc')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('eigval_expo')), 39, 'decimal') line = fillstr(line,
str(p.get_parameter('eigval_expo')), 39, 'decimal')
content += line content += line
line = create_line("N_MAX,N_ITER,N_TABLE,SHIFT") line = create_line("N_MAX,N_ITER,N_TABLE,SHIFT")
line = fillstr(line, str(p.get_parameter('eigval_nmax')), 9, 'left') line = fillstr(line, str(p.get_parameter('eigval_nmax')), 9, 'left')
@ -673,7 +716,8 @@ class SpecIO(object):
thfwd_arr = np.ones((nat)) thfwd_arr = np.ones((nat))
path_filtering = p.extra_parameters['calculation'].get_parameter( path_filtering = p.extra_parameters['calculation'].get_parameter(
'path_filtering').value 'path_filtering').value
if path_filtering != None and 'backward_scattering' in path_filtering: if (path_filtering is not None and
'backward_scattering' in path_filtering):
ibwd_arr = np.ones((nat), dtype=np.int) ibwd_arr = np.ones((nat), dtype=np.int)
else: else:
ibwd_arr = np.zeros((nat), dtype=np.int) ibwd_arr = np.zeros((nat), dtype=np.int)
@ -694,7 +738,8 @@ class SpecIO(object):
line = create_line("IPW,NCUT,PCTINT,IPP") line = create_line("IPW,NCUT,PCTINT,IPP")
line = fillstr(line, str(p.get_parameter('calc_ipw')), 9, 'left') line = fillstr(line, str(p.get_parameter('calc_ipw')), 9, 'left')
line = fillstr(line, str(p.get_parameter('calc_ncut')), 19, 'left') line = fillstr(line, str(p.get_parameter('calc_ncut')), 19, 'left')
line = fillstr(line, str(p.get_parameter('calc_pctint')), 29, 'decimal') line = fillstr(line,
str(p.get_parameter('calc_pctint')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('calc_ipp')), 39, 'left') line = fillstr(line, str(p.get_parameter('calc_ipp')), 39, 'left')
content += line content += line
line = create_line("ILENGTH,RLENGTH,UNLENGTH") line = create_line("ILENGTH,RLENGTH,UNLENGTH")
@ -713,7 +758,8 @@ class SpecIO(object):
idcm = format(2, 'd') idcm = format(2, 'd')
temp = format(0., '.2f') temp = format(0., '.2f')
td = format(500., '.2f') td = format(500., '.2f')
LOGGER.warning('Vibrational damping is disabled for this calculation.') LOGGER.warning('Vibrational damping is disabled for this '
'calculation.')
else: else:
idwsph = str(p.get_parameter('calc_idwsph')) idwsph = str(p.get_parameter('calc_idwsph'))
idcm = str(p.get_parameter('calc_idcm')) idcm = str(p.get_parameter('calc_idcm'))
@ -757,7 +803,8 @@ class SpecIO(object):
content += title(' ' * 17 + 'INPUT FILES (PHD, EXAFS, LEED, AED, ' content += title(' ' * 17 + 'INPUT FILES (PHD, EXAFS, LEED, AED, '
'APECS) :', center=False) 'APECS) :', center=False)
content += rule(tabs=(), symbol='-') content += rule(tabs=(), symbol='-')
content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', content += title(
' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE',
center=False) center=False)
content += rule(tabs=(5, 23, 7, 10)) content += rule(tabs=(5, 23, 7, 10))
line = create_line("DATA FILE,UNIT") line = create_line("DATA FILE,UNIT")
@ -791,7 +838,8 @@ class SpecIO(object):
center=False) center=False)
content += title(' ' * 28 + '(AUGER ELECTRON)', center=False) content += title(' ' * 28 + '(AUGER ELECTRON)', center=False)
content += rule(tabs=(), symbol='-') content += rule(tabs=(), symbol='-')
content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', content += title(
' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE',
center=False) center=False)
content += rule(tabs=(5, 23, 7, 10)) content += rule(tabs=(5, 23, 7, 10))
line = create_line("PHASE SHIFTS/TL FILE,UNIT") line = create_line("PHASE SHIFTS/TL FILE,UNIT")
@ -810,7 +858,8 @@ class SpecIO(object):
content += title(' ' * 29 + 'OUTPUT FILES :', center=False) content += title(' ' * 29 + 'OUTPUT FILES :', center=False)
content += rule(tabs=(), symbol='-') content += rule(tabs=(), symbol='-')
content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', content += title(
' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE',
center=False) center=False)
content += rule(tabs=(5, 23, 7, 10)) content += rule(tabs=(5, 23, 7, 10))
line = create_line("CONTROL FILE,UNIT") line = create_line("CONTROL FILE,UNIT")
@ -826,7 +875,8 @@ class SpecIO(object):
line = fillstr(line, str(p.get_parameter('output_unit11')), 39, 'left') line = fillstr(line, str(p.get_parameter('output_unit11')), 39, 'left')
content += line content += line
line = create_line("AUGMENTED CLUSTER FILE,UNIT") line = create_line("AUGMENTED CLUSTER FILE,UNIT")
line = fillstr(line, str(p.get_parameter('output_augclus')), 9, 'right') line = fillstr(line,
str(p.get_parameter('output_augclus')), 9, 'right')
line = fillstr(line, str(p.get_parameter('output_unit12')), 39, 'left') line = fillstr(line, str(p.get_parameter('output_unit12')), 39, 'left')
content += line content += line
content += rule(tabs=(5, 23, 7, 10)) content += rule(tabs=(5, 23, 7, 10))
@ -874,15 +924,13 @@ class SpecIO(object):
# backup the content in memory # backup the content in memory
old_content = content old_content = content
""" # for key in ('NATP_M', 'NATCLU_M', 'NE_M', 'NEMET_M', 'LI_M', 'NL_M',
for key in ('NATP_M', 'NATCLU_M', 'NE_M', 'NEMET_M', 'LI_M', 'NL_M', # 'NO_ST_M'):
'NO_ST_M'): # required = requirements[key]
required = requirements[key] # limit = self.malloc_parameters.get_parameter(key).value
limit = self.malloc_parameters.get_parameter(key).value # value = required if required > limit else limit
value = required if required > limit else limit # content = re.sub(r'({:s}\s*=\s*)\d+'.format(key),
content = re.sub(r'({:s}\s*=\s*)\d+'.format(key), # r'\g<1>{:d}'.format(value), content)
r'\g<1>{:d}'.format(value), content)
"""
for key in ('NAT_EQ_M', 'N_CL_N_M', 'NDIF_M', 'NSO_M', 'NTEMP_M', for key in ('NAT_EQ_M', 'N_CL_N_M', 'NDIF_M', 'NSO_M', 'NTEMP_M',
'NODES_EX_M', 'NSPIN_M', 'NTH_M', 'NPH_M', 'NDIM_M', 'NODES_EX_M', 'NSPIN_M', 'NTH_M', 'NPH_M', 'NDIM_M',
@ -895,7 +943,6 @@ class SpecIO(object):
content = re.sub(r'({:s}\s*=\s*)\d+'.format(key), content = re.sub(r'({:s}\s*=\s*)\d+'.format(key),
r'\g<1>{:d}'.format(value), content) r'\g<1>{:d}'.format(value), content)
modified = False modified = False
if content != old_content: if content != old_content:
with open(filename, 'w') as fd: with open(filename, 'w') as fd:
@ -945,7 +992,7 @@ class SpecIO(object):
data = np.loadtxt(filename, skiprows=skip, unpack=True) data = np.loadtxt(filename, skiprows=skip, unpack=True)
if len(data.shape) <= 1: if len(data.shape) <= 1:
data = data.reshape((2,-1)) data = data.reshape((2, -1))
return data return data
def load_facdif(self, filename='facdif1.dat'): def load_facdif(self, filename='facdif1.dat'):
@ -956,8 +1003,4 @@ class SpecIO(object):
pat = re.compile(r'ORDER\s+(\d+)\s+TOTAL NUMBER OF PATHS\s+:\s+(\d+)') pat = re.compile(r'ORDER\s+(\d+)\s+TOTAL NUMBER OF PATHS\s+:\s+(\d+)')
with open(filename, 'r') as fd: with open(filename, 'r') as fd:
content = fd.read() content = fd.read()
#return pat.findall(content.replace('\n', '__cr__'))
return pat.findall(content) return pat.findall(content)

151
src/msspec/calculator.py
View File

@ -1,5 +1,25 @@
# coding: utf-8 #!/usr/bin/env python
# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a tw=80: #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/calculator.py
# Last modified: ven. 10 avril 2020 17:19:24
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
""" """
Module calculator Module calculator
================= =================
@ -21,54 +41,58 @@ For more information on MsSpec, follow this
""" """
import inspect
import os import os
import re
import shutil import shutil
import sys import sys
import re
import inspect
from subprocess import Popen, PIPE
from shutil import copyfile, rmtree
from datetime import datetime
import time import time
from io import StringIO
from collections import OrderedDict from collections import OrderedDict
from datetime import datetime
from io import StringIO
from shutil import copyfile
from shutil import rmtree
from subprocess import PIPE
from subprocess import Popen
from ase.calculators.calculator import Calculator
import ase.data
import ase.atom import ase.atom
import ase.atoms import ase.atoms
import ase.data
import numpy as np import numpy as np
from ase.calculators.calculator import Calculator
from terminaltables.ascii_table import AsciiTable
from msspec import iodata from msspec import iodata
from msspec.data import electron_be from msspec.calcio import PhagenIO
from msspec.calcio import SpecIO
from msspec.config import Config from msspec.config import Config
from msspec.misc import (UREG, LOGGER, get_call_info, from msspec.data import electron_be
get_level_from_electron_configuration, from msspec.misc import get_call_info
XRaySource, set_log_output, log_process_output) from msspec.misc import get_level_from_electron_configuration
from msspec.utils import get_atom_index from msspec.misc import log_process_output
from msspec.misc import LOGGER
from msspec.parameters import (PhagenParameters, PhagenMallocParameters, from msspec.misc import set_log_output
SpecParameters, SpecMallocParameters, from msspec.misc import UREG
GlobalParameters, MuffintinParameters, from msspec.misc import XRaySource
TMatrixParameters, SourceParameters, from msspec.parameters import CalculationParameters
DetectorParameters, ScanParameters, from msspec.parameters import DetectorParameters
CalculationParameters, from msspec.parameters import EIGParameters
PEDParameters, EIGParameters) from msspec.parameters import GlobalParameters
from msspec.calcio import PhagenIO, SpecIO from msspec.parameters import MuffintinParameters
from msspec.parameters import PEDParameters
from msspec.parameters import PhagenMallocParameters
from msspec.parameters import PhagenParameters
from msspec.parameters import ScanParameters
from msspec.parameters import SourceParameters
from msspec.parameters import SpecMallocParameters
from msspec.parameters import SpecParameters
from msspec.parameters import TMatrixParameters
from msspec.phagen.fortran.libphagen import main as do_phagen 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_mi
from msspec.spec.fortran import _eig_pw from msspec.spec.fortran import _eig_pw
from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
from msspec.spec.fortran import _phd_se_noso_nosp_nosym
from terminaltables.ascii_table import AsciiTable from msspec.utils import get_atom_index
class _MSCALCULATOR(Calculator): class _MSCALCULATOR(Calculator):
@ -277,34 +301,33 @@ class _MSCALCULATOR(Calculator):
def run_phagen(self): def run_phagen(self):
#include_fname = os.path.join(self.tmp_folder, 'src/msxas3.inc')
input_fname = os.path.join(self.tmp_folder, 'input/input.ms') input_fname = os.path.join(self.tmp_folder, 'input/input.ms')
modified = self.phagenio.write_input_file(filename=input_fname)
#mod0 = self.phagenio.write_include_file(filename=include_fname) # Write the external potential file
mod1 = self.phagenio.write_input_file(filename=input_fname) pot = self.tmatrix_parameters.potential
if pot != 'muffin_tin':
inpot_fname = os.path.join(self.tmp_folder, 'output/fort.2')
mflag = pot.write(inpot_fname, self.phagenio.prototypical_atoms)
modified = modified or mflag
# Make a copy of the input potential file
dest = os.path.join(self.tmp_folder, 'input/extpot.dat')
shutil.copy(inpot_fname, dest)
self.modified = self.modified or mod1 #or mod0 or mod1 # Run if input files changed
if modified:
#if self.modified: # Run phagen
if mod1:
# run phagen
#self._make('tmatrix')
os.chdir(os.path.join(self.tmp_folder, 'output')) 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() do_phagen()
# rename some output files to be more explicit # Copy 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('clus/clus.out', 'cluster.clu') shutil.copy('clus/clus.out', 'cluster.clu')
shutil.copy('fort.35', 'tmatrix.tl') shutil.copy('fort.35', 'tmatrix.tl')
shutil.copy('fort.55', 'tmatrix.rad') shutil.copy('fort.55', 'tmatrix.rad')
# Load data
self.phagenio.load_tl_file('tmatrix.tl')
self.phagenio.load_cluster_file('cluster.clu')
self.phagenio.load_prototypical_atoms('div/molinpot3.out')
# Change back to the init_folder
os.chdir(self.init_folder) os.chdir(self.init_folder)
@ -425,15 +448,19 @@ class _MSCALCULATOR(Calculator):
LOGGER.info("Getting the TMatrix...") LOGGER.info("Getting the TMatrix...")
LOGGER.debug(get_call_info(inspect.currentframe())) LOGGER.debug(get_call_info(inspect.currentframe()))
# If using an external potential, Phagen should be run twice,
# the first time with the internal MT potential.
pot = self.tmatrix_parameters.potential
if pot != 'muffin_tin':
LOGGER.info("Phagen is running first with MuffinTin potential...")
self.tmatrix_parameters.potential = "muffin_tin"
self.tmatrix_parameters.exchange_correlation = "hedin_lundqvist_complex"
self.run_phagen()
self.tmatrix_parameters.potential = pot
self.run_phagen() self.run_phagen()
filename = os.path.join(self.tmp_folder, 'output/tmatrix.tl') tl = self.phagenio.tl
tl = self.phagenio.load_tl_file(filename)
filename = os.path.join(self.tmp_folder, 'output/cluster.clu')
self.phagenio.load_cluster_file(filename)
tl_threshold = self.tmatrix_parameters.get_parameter('tl_threshold') tl_threshold = self.tmatrix_parameters.get_parameter('tl_threshold')
if tl_threshold.value != None: if tl_threshold.value != None:
LOGGER.debug(" applying tl_threshold to %s...", LOGGER.debug(" applying tl_threshold to %s...",

118
src/msspec/config.py
View File

@ -1,118 +0,0 @@
# coding: utf-8
"""
Module config
=============
"""
from configparser import NoSectionError
from configparser import SafeConfigParser as ConfigParser
import os
import sys
import fnmatch
from msspec.misc import LOGGER
class NoConfigFile(Exception): pass
class Config(object):
def __init__(self, **kwargs):
self.fname = os.path.join(os.environ['HOME'], '.config/msspec/pymsspec.cfg')
self.version = sys.modules['msspec'].__version__
self.config = ConfigParser()
self.defaults = {'path': 'None'}
self.defaults.update(kwargs)
# try to load the config file, create one with defaults if none is found
try:
self.options = self.read()
self.options.update(kwargs)
self.set(**self.options)
self.write()
except NoConfigFile:
# write a file with default options
self.write(defaults=True)
except NoSectionError:
# the file exists but has no options for this version of pymsspec
self.config.add_section(self.version)
self.set(**self.defaults)
self.write()
def read(self):
fp = self.config.read(self.fname)
if len(fp) == 0:
raise NoConfigFile
self.version = self.config.get('general', 'version')
return dict(self.config.items(self.version))
def set(self, **kwargs):
if not(self.config.has_section(self.version)):
self.config.add_section(self.version)
for k, v in list(kwargs.items()):
self.config.set(self.version, k, v)
def get(self, key):
return self.config.get(self.version, key)
def choose_msspec_folder(self, *folders):
print("Several folders containing the appropriate version of MsSpec were found.")
print("Please choose one tu use:")
s = ""
for i, f in enumerate(folders):
s += '{:d}) {:s}\n'.format(i, f)
print(s)
return(folders[int(input('Your choice: '))])
def find_msspec_folders(self):
version = sys.modules['msspec'].__version__
folders = []
i = 0
prompt = 'Please wait while scanning the filesystem (%3d found) '
sys.stdout.write(prompt % len(folders))
sys.stdout.write('\033[s')
for root, dirnames, filenames in os.walk('/home/stricot'):
sys.stdout.write('\033[u\033[k')
sys.stdout.write('%d folders scanned' % i)
i += 1
for fn in fnmatch.filter(filenames, 'VERSION'):
with open(os.path.join(root, fn), 'r') as fd:
try:
line = fd.readline()
if line.strip() == version:
folders.append(root)
sys.stdout.write('\033[2000D\033[k')
sys.stdout.write(prompt % len(folders))
except:
pass
sys.stdout.write('\033[u\033[k\n')
print('Done.')
return(folders)
def write(self, defaults=False):
if defaults:
self.set(**self.defaults)
with open(self.fname, 'w') as fd:
self.config.write(fd)
LOGGER.info("{} file written".format(self.fname))
def set_mode(self, mode="pymsspec"):
if not(self.config.has_section("general")):
self.config.add_section("general")
self.config.set("general", "mode", str(mode))
def get_mode(self):
return self.config.get("general", "mode")
def set_version(self, version=""):
if not(self.config.has_section("general")):
self.config.add_section("general")
self.config.set("general", "version", version)
def remove_version(self, version):
pass

26
src/msspec/create_tests_results.py
View File

@ -1,2 +1,28 @@
#!/usr/bin/env python
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/create_tests_results.py
# Last modified: ven. 10 avril 2020 17:29:16
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
from msspec.tests import create_tests_results from msspec.tests import create_tests_results
create_tests_results() create_tests_results()

51
src/msspec/iodata.py
View File

@ -1,4 +1,25 @@
# coding: utf-8 #!/usr/bin/env python
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/iodata.py
# Last modified: ven. 10 avril 2020 17:23:11
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
""" """
Module iodata Module iodata
============= =============
@ -47,29 +68,25 @@ Here is an example of how to store values in a Data object:
import os import os
import numpy as np import sys
import h5py from distutils.version import LooseVersion
from lxml import etree from distutils.version import StrictVersion
import msspec
from msspec.misc import LOGGER
import ase.io
from io import StringIO from io import StringIO
import wx import ase.io
import h5py
import numpy as np
import wx.grid import wx.grid
from lxml import etree
from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas
from matplotlib.backends.backend_wxagg import NavigationToolbar2WxAgg from matplotlib.backends.backend_wxagg import NavigationToolbar2WxAgg
#from matplotlib.backends.backend_wxcairo import FigureCanvasWxCairo as FigureCanvas
#from matplotlib.backends.backend_wxcairo import NavigationToolbar2WxCairo
from matplotlib.figure import Figure from matplotlib.figure import Figure
from terminaltables import AsciiTable from terminaltables import AsciiTable
from distutils.version import StrictVersion, LooseVersion
import sys import msspec
#sys.path.append('../../MsSpecGui/msspecgui/msspec/gui')
from msspec.msspecgui.msspec.gui.clusterviewer import ClusterViewer from msspec.msspecgui.msspec.gui.clusterviewer import ClusterViewer
from msspec.misc import LOGGER
def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1): def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1):
# mix the values of existing theta and new theta and return the # mix the values of existing theta and new theta and return the
@ -1186,7 +1203,3 @@ if __name__ == "__main__":
import sys import sys
data = Data.load(sys.argv[1]) data = Data.load(sys.argv[1])
data.view() data.view()

35
src/msspec/misc.py
View File

@ -1,10 +1,34 @@
# coding: utf-8 #!/usr/bin/env python
# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : # #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/misc.py
# Last modified: ven. 10 avril 2020 17:30:42
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
""" """
Module misc Module misc
=========== ===========
""" """
import logging import logging
from pint import UnitRegistry from pint import UnitRegistry
import numpy as np import numpy as np
@ -12,6 +36,7 @@ import inspect
import re import re
import os import os
class XRaySource(object): class XRaySource(object):
MG_KALPHA = 1253.6 MG_KALPHA = 1253.6
AL_KALPHA = 1486.6 AL_KALPHA = 1486.6
@ -19,10 +44,10 @@ class XRaySource(object):
pass pass
UREG = UnitRegistry() UREG = UnitRegistry()
UREG.define('rydberg = c * h * rydberg_constant = Ry') #UREG.define('rydberg = c * h * rydberg_constant = Ry')
UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0') #UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0')
logging.basicConfig(level=logging.INFO) #logging.basicConfig(level=logging.INFO)
LOGGER = logging.getLogger('msspec') LOGGER = logging.getLogger('msspec')
np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5) np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5)

283
src/msspec/parameters.py
View File

@ -1,16 +1,47 @@
# coding: utf-8 #!/usr/bin/env python
# vim: set et ts=4 sw=4 nu fdm=indent: #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/parameters.py
# Last modified: ven. 10 avril 2020 17:31:50
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
""" """
Module parameters Module parameters
================= =================
""" """
import textwrap
import numpy as np
import re import re
from terminaltables import AsciiTable import textwrap
from msspec.misc import LOGGER, UREG, XRaySource, get_level_from_electron_configuration
import ase import ase
import numpy as np
from terminaltables import AsciiTable
from msspec.misc import get_level_from_electron_configuration
from msspec.misc import LOGGER
from msspec.misc import UREG
from msspec.misc import XRaySource
from msspec.utils import ForeignPotential
from msspec.utils import SPRKKRPotential
class Parameter(object): class Parameter(object):
def __init__(self, name, types=None, limits=(None, None), def __init__(self, name, types=None, limits=(None, None),
@ -31,7 +62,7 @@ class Parameter(object):
self.fmt = fmt self.fmt = fmt
self.binding = binding self.binding = binding
self.private = private self.private = private
self.docstring = doc self.docstring = textwrap.dedent(doc)
self._value = None self._value = None
self.value = default self.value = default
@ -59,9 +90,7 @@ class Parameter(object):
return value.to(self.unit).magnitude return value.to(self.unit).magnitude
return value return value
def check(self, value): def assert_message(self, msg, *args):
def assert_message(msg, *args):
s = '\'{}\': {}'.format(self.name, msg) s = '\'{}\': {}'.format(self.name, msg)
s = s.format(*args) s = s.format(*args)
@ -94,6 +123,8 @@ class Parameter(object):
return s return s
def check(self, value):
# convert if a unit was given # convert if a unit was given
_value = self.convert(value) _value = self.convert(value)
@ -106,7 +137,7 @@ class Parameter(object):
assert bool in self.allowed_types assert bool in self.allowed_types
assert isinstance(_value, self.allowed_types) assert isinstance(_value, self.allowed_types)
except AssertionError: except AssertionError:
raise AssertionError(assert_message( raise AssertionError(self.assert_message(
'Type {} is not an allowed type for this option ' 'Type {} is not an allowed type for this option '
'({} allowed)', '({} allowed)',
str(type(_value)), str(self.allowed_types))) str(type(_value)), str(self.allowed_types)))
@ -120,22 +151,22 @@ class Parameter(object):
# for i, val in enumerate(_values): # for i, val in enumerate(_values):
for val in np.array(_value).flatten(): for val in np.array(_value).flatten():
if self.low_limit != None: if self.low_limit != None:
assert val >= self.low_limit, assert_message( assert val >= self.low_limit, self.assert_message(
'Value must be >= {:s}', 'Value must be >= {:s}',
str(self.low_limit)) str(self.low_limit))
if self.high_limit != None: if self.high_limit != None:
assert val <= self.high_limit, assert_message( assert val <= self.high_limit, self.assert_message(
'Value must be <= {:s}', 'Value must be <= {:s}',
str(self.high_limit)) str(self.high_limit))
if self.allowed_values != None and isinstance(val, tuple(type(x) for x in self.allowed_values)): if self.allowed_values != None: # and isinstance(val, tuple(type(x) for x in self.allowed_values)):
assert val in self.allowed_values, assert_message( assert val in self.allowed_values, self.assert_message(
'This value is not allowed. Please choose between ' 'This value is not allowed. Please choose between '
'one of {:s}', 'one of {:s}',
str(self.allowed_values)) str(self.allowed_values))
if self.pattern != None: if self.pattern != None:
p = re.compile(self.pattern) p = re.compile(self.pattern)
m = p.match(val) m = p.match(val)
assert m != None, assert_message( assert m != None, self.assert_message(
'Wrong string format') 'Wrong string format')
# _value[i] = val # _value[i] = val
@ -157,8 +188,8 @@ class Parameter(object):
new_value = None new_value = None
try: try:
new_value = self.binding(self) new_value = self.binding(self)
except AttributeError: except AttributeError as err:
pass LOGGER.warning(err)
if new_value is not None: if new_value is not None:
self._value = new_value self._value = new_value
# LOGGER.debug('{:s} = {:s}'.format(self.name, str(self.value))) # LOGGER.debug('{:s} = {:s}'.format(self.name, str(self.value)))
@ -716,7 +747,7 @@ class GlobalParameters(BaseParameters):
parameters = ( parameters = (
Parameter('spectroscopy', types=str, allowed_values=( Parameter('spectroscopy', types=str, allowed_values=(
'PED', 'AED', 'LEED', 'EXAFS', 'EIG'), default='PED', 'PED', 'AED', 'LEED', 'EXAFS', 'EIG'), default='PED',
doc=textwrap.dedent(""" doc="""
There are 4 choices for the spectroscopy option: There are 4 choices for the spectroscopy option:
- '**PED**', for Photo Electron Diffraction - '**PED**', for Photo Electron Diffraction
@ -729,12 +760,12 @@ class GlobalParameters(BaseParameters):
matrix. matrix.
The value is case insensitive. The value is case insensitive.
""")), """),
Parameter('algorithm', types=str, allowed_values=('expansion', Parameter('algorithm', types=str, allowed_values=('expansion',
'inversion', 'inversion',
'correlation', 'correlation',
'power'), 'power'),
default='expansion', doc=textwrap.dedent(""" default='expansion', doc="""
You can choose the algorithm used for the computation of the scattering path operator. You can choose the algorithm used for the computation of the scattering path operator.
- '**inversion**', for the classical matrix inversion method - '**inversion**', for the classical matrix inversion method
@ -745,10 +776,10 @@ class GlobalParameters(BaseParameters):
The series expansion algorithm is well suited for high energy since the number of terms The series expansion algorithm is well suited for high energy since the number of terms
required decreases as the energy increases. The exact solution is obtained by the matrix inversion required decreases as the energy increases. The exact solution is obtained by the matrix inversion
method but should be preferably used for lower energy. method but should be preferably used for lower energy.
""")), """),
Parameter('polarization', types=(type(None), str), Parameter('polarization', types=(type(None), str),
allowed_values=(None, 'linear_qOz', 'linear_xOy', allowed_values=(None, 'linear_qOz', 'linear_xOy',
'circular'), default=None, doc=textwrap.dedent(""" 'circular'), default=None, doc="""
Specify the polarization of the incident light. Specify the polarization of the incident light.
- **None**, for unpolarized light - **None**, for unpolarized light
@ -756,16 +787,16 @@ class GlobalParameters(BaseParameters):
- '**linear_xOy**' for a polarization vector in the :math:`(x0y)` plane - '**linear_xOy**' for a polarization vector in the :math:`(x0y)` plane
- '**circular**' for circular dichroism - '**circular**' for circular dichroism
""")), """),
Parameter('dichroism', types=(type(None), str), Parameter('dichroism', types=(type(None), str),
allowed_values=(None, 'natural', 'sum_over_spin', allowed_values=(None, 'natural', 'sum_over_spin',
'spin_resolved'), default=None, doc=textwrap.dedent(""" 'spin_resolved'), default=None, doc="""
Used to perform dichroic calculations. The default (None) is to disable this. Used to perform dichroic calculations. The default (None) is to disable this.
""")), """),
Parameter('spinpol', types=bool, default=False, doc=textwrap.dedent(""" Parameter('spinpol', types=bool, default=False, doc="""
Enable or disbale spin-resolved calculations. Enable or disbale spin-resolved calculations.
""")), """),
Parameter('folder', types=str, default='./calc', doc=textwrap.dedent(""" Parameter('folder', types=str, default='./calc', doc="""
This parameter is the path to the temporary folder used for the calculations. If you do not change this This parameter is the path to the temporary folder used for the calculations. If you do not change this
parameter between calculations, all the content will be overridden. This is usually not a problem, since the parameter between calculations, all the content will be overridden. This is usually not a problem, since the
whole bunch of data created during a computation is not meant to be saved. But you may want to anyway by whole bunch of data created during a computation is not meant to be saved. But you may want to anyway by
@ -781,7 +812,7 @@ class GlobalParameters(BaseParameters):
calc.shutdown() # the folder is removed calc.shutdown() # the folder is removed
""")) """)
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -816,11 +847,11 @@ class GlobalParameters(BaseParameters):
class MuffintinParameters(BaseParameters): class MuffintinParameters(BaseParameters):
def __init__(self, phagen_parameters, spec_parameters): def __init__(self, phagen_parameters, spec_parameters):
parameters = ( parameters = (
Parameter('charge_relaxation', types=bool, default=True, doc=textwrap.dedent(""" Parameter('charge_relaxation', types=bool, default=True, doc="""
Used to specify wether the charge density of the photoabsorbing atom, which is used Used to specify wether the charge density of the photoabsorbing atom, which is used
to construct the potential, is allowaed to relax around the core hole or not. to construct the potential, is allowaed to relax around the core hole or not.
""")), """),
Parameter('ionicity', types=dict, default={}, doc=textwrap.dedent(""" Parameter('ionicity', types=dict, default={}, doc="""
A dictionary to specify the ionicity of each kind of atoms. If empty, the atoms are considered to be A dictionary to specify the ionicity of each kind of atoms. If empty, the atoms are considered to be
neutral. Otherwise, each key must be a chemical symbol and the corresponding value should be the fraction neutral. Otherwise, each key must be a chemical symbol and the corresponding value should be the fraction
of electrons added (negative) or substracted (positive) with respect to neutrality. of electrons added (negative) or substracted (positive) with respect to neutrality.
@ -831,28 +862,28 @@ class MuffintinParameters(BaseParameters):
means that 0.15 electrons have been substracted from La, likewise from Fe. Neutrality implies that 0.3 means that 0.15 electrons have been substracted from La, likewise from Fe. Neutrality implies that 0.3
electrons have to be added to oxygen atoms. electrons have to be added to oxygen atoms.
""")), """),
Parameter('relativistic_mode', types=str, Parameter('relativistic_mode', types=str,
allowed_values=('non_relativistic', 'scalar_relativistic', allowed_values=('non_relativistic', 'scalar_relativistic',
'spin_orbit_resolved'), 'spin_orbit_resolved'),
default='non_relativistic', doc=textwrap.dedent(""" default='non_relativistic', doc="""
To tell whether to use the scalar relativistic approximation or not. To tell whether to use the scalar relativistic approximation or not.
""")), """),
Parameter('radius_overlapping', types=float, limits=(0., 1.), Parameter('radius_overlapping', types=float, limits=(0., 1.),
default=0., doc=textwrap.dedent(""" default=0., doc="""
to allow atomic spheres to overlapp with each other. The value is a percentage, 1. means 100%. to allow atomic spheres to overlapp with each other. The value is a percentage, 1. means 100%.
""")), """),
Parameter('interstitial_potential', types=(int, float), Parameter('interstitial_potential', types=(int, float),
unit=UREG.eV, default=0., doc=textwrap.dedent(""" unit=UREG.eV, default=0., doc="""
The average interstitial potential (or inner potential) expressed in eV. It is used to compute The average interstitial potential (or inner potential) expressed in eV. It is used to compute
the refraction at the surface. the refraction at the surface.
""")), """),
Parameter('hydrogen_radius', types=(int, float), default=0.9, Parameter('hydrogen_radius', types=(int, float), default=0.9,
limits=(0., None), unit=UREG.angstroms, doc=textwrap.dedent(""" limits=(0., None), unit=UREG.angstroms, doc="""
The program can have difficulties to find the radius of the hydrogen atom (small atom). You can The program can have difficulties to find the radius of the hydrogen atom (small atom). You can
specify here a value for the radius. If you set it to 'None', the calculation of the muffin-tin specify here a value for the radius. If you set it to 'None', the calculation of the muffin-tin
radius of H atoms will be left to the program. radius of H atoms will be left to the program.
""")) """)
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -884,28 +915,28 @@ class MuffintinParameters(BaseParameters):
class TMatrixParameters(BaseParameters): class TMatrixParameters(BaseParameters):
def __init__(self, phagen_parameters): def __init__(self, phagen_parameters):
parameters = ( parameters = (
Parameter('potential', types=str, Parameter('potential', types=(str, ForeignPotential),
allowed_values=('muffin_tin', 'lmto', 'spkkr', 'msf'), default='muffin_tin',
default='muffin_tin', doc=textwrap.dedent(""" doc="""
This option allows to choose which kind of potential is used to compute the T-Matrix. For now, This option allows to choose which kind of potential is
only the internal *Muffin-Tin* potential is supported. 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, Parameter('exchange_correlation', types=str,
allowed_values=('hedin_lundqvist_real', allowed_values=('hedin_lundqvist_real',
'hedin_lundqvist_complex', 'hedin_lundqvist_complex',
'x_alpha_real', 'x_alpha_real',
'dirac_hara_real', 'dirac_hara_complex'), 'dirac_hara_real', 'dirac_hara_complex'),
default='hedin_lundqvist_complex', doc=textwrap.dedent(""" default='hedin_lundqvist_complex', doc="""
Set the type of exchange and correlation potential that will be used. Set the type of exchange and correlation potential that will be used.
""")), """),
Parameter('imaginery_part', types=(int, float), default=0., doc=textwrap.dedent(""" Parameter('imaginery_part', types=(int, float), default=0., doc="""
This value is added to complex potentials to account for core-hole lifetime and experimental resolution This value is added to complex potentials to account for core-hole lifetime and experimental resolution
broadening effects. broadening effects.
""")), """),
Parameter('lmax_mode', types=str, Parameter('lmax_mode', types=str,
allowed_values=('max_ke', 'true_ke', 'imposed'), allowed_values=('max_ke', 'true_ke', 'imposed'),
default='true_ke', doc=textwrap.dedent(""" default='true_ke', doc="""
This allows to control the number of basis functions used to expand the wave function on each This allows to control the number of basis functions used to expand the wave function on each
atom. It can be: atom. It can be:
@ -919,11 +950,11 @@ class TMatrixParameters(BaseParameters):
:math:`l_{max} = kr_{at} + 1` where :math:`k=E^{1/2}_k` with :math:`E_k` being the kinetic :math:`l_{max} = kr_{at} + 1` where :math:`k=E^{1/2}_k` with :math:`E_k` being the kinetic
energy. In this case :math:`l_{max}` depends both on the energy and the atom number. energy. In this case :math:`l_{max}` depends both on the energy and the atom number.
""")), """),
Parameter('lmaxt', types=int, limits=(1, None), default=19, doc=textwrap.dedent(""" Parameter('lmaxt', types=int, limits=(1, None), default=19, doc="""
The value of :math:`l_{max}` if *lmax_mode = 'imposed'* The value of :math:`l_{max}` if *lmax_mode = 'imposed'*
""")), """),
Parameter('tl_threshold', types=(type(None), float), default=None, doc=textwrap.dedent(""" Parameter('tl_threshold', types=(type(None), float), default=None, doc="""
This option allows to control the number of basis function by defining a threshold value for the *tl*. This option allows to control the number of basis function by defining a threshold value for the *tl*.
For example:: For example::
@ -934,8 +965,8 @@ class TMatrixParameters(BaseParameters):
.. note:: .. note::
This option is compatible with any modes of the *lmax_mode* option. This option is compatible with any modes of the *lmax_mode* option.
""")), """),
Parameter('max_tl', types=(type(None), dict), default=None, doc=textwrap.dedent(""" Parameter('max_tl', types=(type(None), dict), default=None, doc="""
*max_tl* is used to sepcify a maximum number of basis functions to use for each kind of atoms. For example, *max_tl* is used to sepcify a maximum number of basis functions to use for each kind of atoms. For example,
in the case of an MgO cluster, you could write:: in the case of an MgO cluster, you could write::
@ -948,7 +979,7 @@ class TMatrixParameters(BaseParameters):
.. note:: .. note::
This option is compatible with any modes of the *lmax_mode* option. This option is compatible with any modes of the *lmax_mode* option.
""")) """)
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -956,12 +987,17 @@ class TMatrixParameters(BaseParameters):
self.freeze() self.freeze()
def bind_potential(self, p): def bind_potential(self, p):
mapping = {'muffin_tin': 'in', 'lmto': 'ex', 'spkkr': 'ex', 'msf': 'ex'} mapping = {SPRKKRPotential: 'spkkr'}
if p.value.lower() in ('muffin_tin'): if isinstance(p.value, str):
assert p.value.lower() in ('muffin_tin'), p.assert_message(
"\"{}\" is not an known potential type", p.value)
self.phagen_parameters.potgen = 'in' self.phagen_parameters.potgen = 'in'
elif p.value.lower() in ('spkkr', 'lmto', 'msf'): elif isinstance(p.value, ForeignPotential):
self.phagen_parameters.potgen = 'ex' self.phagen_parameters.potgen = 'ex'
self.phagen_parameters.potype = p.value.lower() potype = mapping[p.value.__class__]
self.phagen_parameters.potype = potype
def bind_exchange_correlation(self, p): def bind_exchange_correlation(self, p):
potential = self.get_parameter('potential').value potential = self.get_parameter('potential').value
@ -975,7 +1011,8 @@ class TMatrixParameters(BaseParameters):
} }
self.phagen_parameters.potype = mapping[p.value] self.phagen_parameters.potype = mapping[p.value]
else: else:
self.phagen_parameters.potype = potential mapping = {SPRKKRPotential: 'spkkr'}
self.phagen_parameters.potype = mapping[potential.value.__class__]
def bind_potential_file(self, p): def bind_potential_file(self, p):
pass pass
@ -1009,7 +1046,7 @@ class SourceParameters(BaseParameters):
def __init__(self, global_parameters=None, phagen_parameters=None, spec_parameters=None): def __init__(self, global_parameters=None, phagen_parameters=None, spec_parameters=None):
parameters = ( parameters = (
Parameter('energy', types=(list, tuple, int, float), Parameter('energy', types=(list, tuple, int, float),
limits=(0, None), unit=UREG.eV, doc=textwrap.dedent(""" limits=(0, None), unit=UREG.eV, doc="""
The photon energy in eV. The photon energy in eV.
Common laboratories X-ray source Mg |kalpha| and Al |kalpha| lines are Common laboratories X-ray source Mg |kalpha| and Al |kalpha| lines are
@ -1026,20 +1063,20 @@ class SourceParameters(BaseParameters):
1253.6 1253.6
"""), """,
default=XRaySource.MG_KALPHA), default=XRaySource.MG_KALPHA),
Parameter('theta', types=(int, float), limits=(-180., 180.), Parameter('theta', types=(int, float), limits=(-180., 180.),
unit=UREG.degree, default=-55., doc=textwrap.dedent(""" unit=UREG.degree, default=-55., doc="""
The polar angle of the photon incidence (in degrees). Please refer to The polar angle of the photon incidence (in degrees). Please refer to
:ref:`this figure <ped_full_picture>` for questions regarding the proper :ref:`this figure <ped_full_picture>` for questions regarding the proper
orientation. orientation.
""")), """),
Parameter('phi', types=(int, float), limits=(-180., 180.), Parameter('phi', types=(int, float), limits=(-180., 180.),
unit=UREG.degree, default=0., doc=textwrap.dedent(""" unit=UREG.degree, default=0., doc="""
The azimuthal angle of the photon incidence (in degrees). Please refer to The azimuthal angle of the photon incidence (in degrees). Please refer to
:ref:`this figure <ped_full_picture>` for questions regarding the proper :ref:`this figure <ped_full_picture>` for questions regarding the proper
orientation. orientation.
""")), """),
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -1091,13 +1128,13 @@ class DetectorParameters(BaseParameters):
parameters = ( parameters = (
Parameter('angular_acceptance', types=(int, float), Parameter('angular_acceptance', types=(int, float),
unit=UREG.degree, limits=(0., None), default=0., unit=UREG.degree, limits=(0., None), default=0.,
doc=textwrap.dedent(""" doc="""
The angular acceptance of the detector in degrees used The angular acceptance of the detector in degrees used
when the *average_sampling* option is enabled below. when the *average_sampling* option is enabled below.
""")), """),
Parameter('average_sampling', types=(type(None), str), Parameter('average_sampling', types=(type(None), str),
allowed_values=(None, 'low', 'medium', 'high'), allowed_values=(None, 'low', 'medium', 'high'),
default=None, doc=textwrap.dedent(""" default=None, doc="""
Used to averaged the signal over directions lying in the Used to averaged the signal over directions lying in the
cone of half-angle *angular_acceptance*. The number of cone of half-angle *angular_acceptance*. The number of
directions to take into account depends on the choosen directions to take into account depends on the choosen
@ -1107,11 +1144,11 @@ class DetectorParameters(BaseParameters):
- '**low**', to average over 5 directions - '**low**', to average over 5 directions
- '**medium**', to average over 13 directions - '**medium**', to average over 13 directions
- '**high**', to average over 49 directions - '**high**', to average over 49 directions
""")), """),
Parameter('rotate', types=bool, default=False, doc=textwrap.dedent(""" Parameter('rotate', types=bool, default=False, doc="""
When False, the sample is rotated when doing a scan (the When False, the sample is rotated when doing a scan (the
usual way). Otherwise, the detector is rotated. usual way). Otherwise, the detector is rotated.
""")) """)
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -1172,7 +1209,7 @@ class ScanParameters(BaseParameters):
default=np.array([0.])), default=np.array([0.])),
Parameter('kinetic_energy', types=(list, tuple, int, float), Parameter('kinetic_energy', types=(list, tuple, int, float),
unit=UREG.eV, limits=(0., None), unit=UREG.eV, limits=(0., None),
default=200., doc=textwrap.dedent(""" default=200., doc="""
if given as a list or tuple: if given as a list or tuple:
* with 2 elements, 10 points of energy will be generated * with 2 elements, 10 points of energy will be generated
with the first element as the starting energy and the with the first element as the starting energy and the
@ -1183,7 +1220,7 @@ class ScanParameters(BaseParameters):
if given as a float or integer, there will be one point if given as a float or integer, there will be one point
for the kinetic energy. for the kinetic energy.
""")), """),
Parameter('ke_array', types=np.ndarray, unit=UREG.eV, Parameter('ke_array', types=np.ndarray, unit=UREG.eV,
default=np.array([200., ]), private=True) default=np.array([200., ]), private=True)
) )
@ -1380,48 +1417,48 @@ class CalculationParameters(BaseParameters):
def __init__(self, global_parameters, phagen_parameters, spec_parameters): def __init__(self, global_parameters, phagen_parameters, spec_parameters):
parameters = ( parameters = (
Parameter('RA_cutoff', types=int, limits=(0, 8), default=1, Parameter('RA_cutoff', types=int, limits=(0, 8), default=1,
doc=textwrap.dedent(""" doc="""
The Rehr-Albers cut-off parameter which controls the degree of The Rehr-Albers cut-off parameter which controls the degree of
sphericity introduced in the description of the basis functions sphericity introduced in the description of the basis functions
used to expand the wave function around each atomic center. used to expand the wave function around each atomic center.
It is only meaningful for the series expansion algorithm. It is only meaningful for the series expansion algorithm.
Its value is limited to 8 but it is rarely necessary to go beyond Its value is limited to 8 but it is rarely necessary to go beyond
2 or 3.""")), 2 or 3."""),
Parameter('scattering_order', types=int, limits=(1, 10), default=3, Parameter('scattering_order', types=int, limits=(1, 10), default=3,
doc=textwrap.dedent(""" doc="""
The scattering order. Only meaningful for the 'expansion' algorithm. The scattering order. Only meaningful for the 'expansion' algorithm.
Its value is limited to 10.""")), 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', 'Pi_1', 'Lowdin'), 'Z_n', 'Pi_1', 'Lowdin'),
types=(type(None), str), default=None, types=(type(None), str), default=None,
doc=textwrap.dedent(""" doc="""
Enable the calculation of the coefficients for the renormalization of Enable the calculation of the coefficients for the renormalization of
the multiple scattering series. the multiple scattering series.
You can choose to renormalize in terms of the :math:`G_n`, the 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:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin
:math:`K^2` matrices""")), :math:`K^2` matrices"""),
Parameter('renormalization_omega', types=(int,float,complex), Parameter('renormalization_omega', types=(int,float,complex),
default=1.+0j, default=1.+0j,
doc=textwrap.dedent(""" doc="""
The :math:`\\omega` coefficient used to initialize the The :math:`\\omega` coefficient used to initialize the
renormalization alogorithm.""")), renormalization alogorithm."""),
Parameter('RA_cutoff_damping', types=int, limits=(0, None), Parameter('RA_cutoff_damping', types=int, limits=(0, None),
default=0, doc=textwrap.dedent(""" default=0, doc="""
The Rehr-Albers truncation order. If > 0, the *RA_cutoff* is The Rehr-Albers truncation order. If > 0, the *RA_cutoff* is
decreased by 1 every *i*\ :sup:`th` scatterer until 0, where decreased by 1 every *i*\ :sup:`th` scatterer until 0, where
*i* = *RA_cutoff_damping*.""")), *i* = *RA_cutoff_damping*."""),
Parameter('spin_flip', types=bool, default=False, Parameter('spin_flip', types=bool, default=False,
doc=textwrap.dedent(""" doc="""
This parameter tells if spin-flip is authorized or not in the This parameter tells if spin-flip is authorized or not in the
scattering process. scattering process.
:Note: :Note:
This option works only if the spin polarization is This option works only if the spin polarization is
enabled in your calculator object (see spinpol_).""")), enabled in your calculator object (see spinpol_)."""),
Parameter('integrals', types=str, allowed_values=('all', Parameter('integrals', types=str, allowed_values=('all',
'diagonal'), 'diagonal'),
default='all', doc=textwrap.dedent(""" default='all', doc="""
This option allows to take into account all four radial integrals This option allows to take into account all four radial integrals
(R++, R+-, R-+ and R--) in the calculation or only the diagonal (R++, R+-, R-+ and R--) in the calculation or only the diagonal
radial integrals (R++ and R--) which are generally much larger. radial integrals (R++ and R--) which are generally much larger.
@ -1431,7 +1468,7 @@ class CalculationParameters(BaseParameters):
This option works only if the spin polarization is This option works only if the spin polarization is
enabled in your calculator object. enabled in your calculator object.
""")), """),
Parameter('path_filtering', types=(type(None), str, tuple, list), Parameter('path_filtering', types=(type(None), str, tuple, list),
allowed_values=(None, allowed_values=(None,
'forward_scattering', 'forward_scattering',
@ -1439,66 +1476,66 @@ class CalculationParameters(BaseParameters):
'distance_cutoff', 'distance_cutoff',
'plane_wave_normal', 'plane_wave_normal',
'plane_wave_spin_averaged'), 'plane_wave_spin_averaged'),
default=None, doc=textwrap.dedent(""" default=None, doc="""
Used to activate some filters. It is possible to specify several Used to activate some filters. It is possible to specify several
of them by grouping them in a tuple or a list. For example:: of them by grouping them in a tuple or a list. For example::
>>> my_filters = ('forward_scattering', 'backward_scattering') >>> my_filters = ('forward_scattering', 'backward_scattering')
>>> calc.calculation_parameters.path_filtering = my_filters >>> calc.calculation_parameters.path_filtering = my_filters
""")), """),
Parameter('off_cone_events', types=int, limits=(0, None), default=1, Parameter('off_cone_events', types=int, limits=(0, None), default=1,
doc=textwrap.dedent(""" doc="""
Used in conjunction with the '*forward_scattering*' filter. Used in conjunction with the '*forward_scattering*' filter.
If the number of scattering processes outside the forward (or If the number of scattering processes outside the forward (or
backward) scattering cone is greater than this number, then the backward) scattering cone is greater than this number, then the
path is rejected and its contribution to the scattering path path is rejected and its contribution to the scattering path
operator won’t be computed. operator won’t be computed.
""")), """),
Parameter('scattering_order_cutoff', types=int, limits=(0, 10), Parameter('scattering_order_cutoff', types=int, limits=(0, 10),
default=2, doc=textwrap.dedent(""" default=2, doc="""
Used in conjunction with the ‘*plane_wave_normal*’ filter. It states Used in conjunction with the ‘*plane_wave_normal*’ filter. It states
to activate the plane wave approximation (which is fast but to activate the plane wave approximation (which is fast but
less accurate) to compute the contribution when the scattering order less accurate) to compute the contribution when the scattering order
is greater than this value.""")), is greater than this value."""),
Parameter('distance', types=(int, float), limits=(0, None), Parameter('distance', types=(int, float), limits=(0, None),
unit=UREG.angstroms, default=10., doc=textwrap.dedent(""" unit=UREG.angstroms, default=10., doc="""
Used with the '*distance_cut_off*' filter. Paths whose length is Used with the '*distance_cut_off*' filter. Paths whose length is
larger than this value are simply rejected.""")), larger than this value are simply rejected."""),
Parameter('vibrational_damping', types=(type(None), str), Parameter('vibrational_damping', types=(type(None), str),
allowed_values=(None, 'debye_waller', 'averaged_tl'), allowed_values=(None, 'debye_waller', 'averaged_tl'),
default='debye_waller', doc=textwrap.dedent(""" default='debye_waller', doc="""
Tells how to compute the effect of atomic vibrations. It can be: Tells how to compute the effect of atomic vibrations. It can be:
- '**debye_waller**' for using the Debye Waller model. - '**debye_waller**' for using the Debye Waller model.
- '**averaged_tl**' to use the more correct averaging over T-matrix elements.""")), - '**averaged_tl**' to use the more correct averaging over T-matrix elements."""),
Parameter('temperature', types=(int, float), limits=(0, None), Parameter('temperature', types=(int, float), limits=(0, None),
unit=UREG.degK, default=293., doc=textwrap.dedent(""" unit=UREG.degK, default=293., doc="""
The temperature of the cluster. Used when *use_debye_model* = True The temperature of the cluster. Used when *use_debye_model* = True
""")), """),
Parameter('debye_temperature', types=(int, float), limits=(0, None), Parameter('debye_temperature', types=(int, float), limits=(0, None),
unit=UREG.degK, default=420., doc=textwrap.dedent(""" unit=UREG.degK, default=420., doc="""
The Debye temperature used for the calculation of the mean square The Debye temperature used for the calculation of the mean square
displacements if *use_debye_model* = True""")), displacements if *use_debye_model* = True"""),
Parameter('use_debye_model', types=bool, default=False, Parameter('use_debye_model', types=bool, default=False,
doc=textwrap.dedent(""" doc="""
No matter the way you compute the effect of atomic vibrations, No matter the way you compute the effect of atomic vibrations,
you need the mean square displacements of atoms. It can be computed you need the mean square displacements of atoms. It can be computed
internally following the Debye model if you set this option to True. internally following the Debye model if you set this option to True.
""")), """),
Parameter('vibration_scaling', types=(int, float), Parameter('vibration_scaling', types=(int, float),
limits=(0., None), default=1.2, doc=textwrap.dedent(""" limits=(0., None), default=1.2, doc="""
Used to simulate the fact that surface atoms vibrate more than Used to simulate the fact that surface atoms vibrate more than
bulk ones. It is a factor applied to the mean square displacements. bulk ones. It is a factor applied to the mean square displacements.
""")), """),
Parameter('basis_functions', types=str, allowed_values=( Parameter('basis_functions', types=str, allowed_values=(
'plane_wave', 'spherical'), default='spherical', private=True), 'plane_wave', 'spherical'), default='spherical', private=True),
Parameter('cutoff_factor', types=(int, float), Parameter('cutoff_factor', types=(int, float),
limits=(1e-4, 999.9999), default=0.01, private=True), limits=(1e-4, 999.9999), default=0.01, private=True),
Parameter('mean_free_path', types=(int, float, str), Parameter('mean_free_path', types=(int, float, str),
default='SeahDench', allowed_values=('mono', 'SeahDench'), default='SeahDench', allowed_values=('mono', 'SeahDench'),
doc=textwrap.dedent(""" doc="""
The electron mean free path value. You can either: The electron mean free path value. You can either:
- Enter a value (in Angströms), in this case any value <=0 will disable the damping - Enter a value (in Angströms), in this case any value <=0 will disable the damping
- Enter the keyword 'mono' to use a formula valid only for monoelemental samples - Enter the keyword 'mono' to use a formula valid only for monoelemental samples
@ -1510,7 +1547,7 @@ class CalculationParameters(BaseParameters):
to a real potential as, when the potential is complex, this damping is taken care to a real potential as, when the potential is complex, this damping is taken care
of by the imaginery part othe potential. of by the imaginery part othe potential.
""")), """),
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
@ -1728,7 +1765,7 @@ class PEDParameters(BaseParameters):
def __init__(self, phagen_parameters, spec_parameters): def __init__(self, phagen_parameters, spec_parameters):
parameters = ( parameters = (
Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$', Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$',
default='1s', doc=textwrap.dedent(""" default='1s', doc="""
The level is the electronic level where the electron comes from. The level is the electronic level where the electron comes from.
It is written: *nlJ* It is written: *nlJ*
where: where:
@ -1742,7 +1779,7 @@ class PEDParameters(BaseParameters):
>>> calc.spectroscopy_parameters.level = '2p3/2' >>> calc.spectroscopy_parameters.level = '2p3/2'
>>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2' >>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2'
""")), """),
Parameter('final_state', types=int, limits=(-1, 2), default=2), Parameter('final_state', types=int, limits=(-1, 2), default=2),
Parameter('spin_orbit', types=(type(None), str), Parameter('spin_orbit', types=(type(None), str),
allowed_values=(None, 'single', 'both'), default=None), allowed_values=(None, 'single', 'both'), default=None),
@ -1779,7 +1816,7 @@ class EIGParameters(BaseParameters):
def __init__(self, phagen_parameters, spec_parameters): def __init__(self, phagen_parameters, spec_parameters):
parameters = ( parameters = (
Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$', Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$',
default='1s', doc=textwrap.dedent(""" default='1s', doc="""
The level is the electronic level where the electron comes from. The level is the electronic level where the electron comes from.
It is written: *nlJ* It is written: *nlJ*
where: where:
@ -1793,21 +1830,21 @@ class EIGParameters(BaseParameters):
>>> calc.spectroscopy_parameters.level = '2p3/2' >>> calc.spectroscopy_parameters.level = '2p3/2'
>>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2' >>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2'
""")), """),
Parameter('final_state', types=int, limits=(-1, 2), default=2), Parameter('final_state', types=int, limits=(-1, 2), default=2),
Parameter('spin_orbit', types=(type(None), str), Parameter('spin_orbit', types=(type(None), str),
allowed_values=(None, 'single', 'both'), default=None), allowed_values=(None, 'single', 'both'), default=None),
Parameter('kernel_matrix_spectrum', types=(bool,), default=False, doc=textwrap.dedent(""" Parameter('kernel_matrix_spectrum', types=(bool,), default=False, doc="""
Whether to output the kernel matrix spectrum for each energy point. Whether to output the kernel matrix spectrum for each energy point.
""")), """),
Parameter('method', types=(str,), default='EPSI', Parameter('method', types=(str,), default='EPSI',
allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG', allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG',
'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER', 'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER',
'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV', 'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV',
'ELEV', 'EULE', 'GBWT', 'VARI', 'ITHE', 'ELEV', 'EULE', 'GBWT', 'VARI', 'ITHE',
'EALG'], 'EALG'],
doc=textwrap.dedent("""The convergence acceleration scheme to be used. doc="""The convergence acceleration scheme to be used.
""")), """),
) )
BaseParameters.__init__(self) BaseParameters.__init__(self)
self.add_parameters(*parameters) self.add_parameters(*parameters)
@ -1832,5 +1869,3 @@ class EIGParameters(BaseParameters):
def bind_kernel_matrix_spectrum(self, p): def bind_kernel_matrix_spectrum(self, p):
value = int(p.value) value = int(p.value)
self.spec_parameters.eigval_ispectrum_ne = value self.spec_parameters.eigval_ispectrum_ne = value

38
src/msspec/tests.py
View File

@ -1,21 +1,41 @@
# coding: utf-8 #!/usr/bin/env python
#
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/tests.py
# Last modified: ven. 10 avril 2020 17:33:28
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
from msspec.calculator import MSSPEC
from msspec.utils import *
from msspec.misc import set_log_level
from ase.build import bulk
import os
import sys
import unittest
from hashlib import md5 from hashlib import md5
from pickle import dumps from pickle import dumps
import unittest from ase.build import bulk
import sys
import os from msspec.calculator import MSSPEC
from msspec.misc import set_log_level
from msspec.utils import *
set_log_level('error') set_log_level('error')
RESULTS_FILENAME = os.path.join(os.path.dirname(__file__), 'results.txt') RESULTS_FILENAME = os.path.join(os.path.dirname(__file__), 'results.txt')
def perform_test(obj, funcname, filename=RESULTS_FILENAME): def perform_test(obj, funcname, filename=RESULTS_FILENAME):
f = getattr(obj, '_' + funcname) f = getattr(obj, '_' + funcname)
output = md5(dumps(f())).hexdigest() output = md5(dumps(f())).hexdigest()

380
src/msspec/utils.py
View File

@ -1,5 +1,26 @@
# -*- encoding: utf-8 -*- #!/usr/bin/env python
# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : # #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/utils.py
# Last modified: ven. 10 avril 2020 15:49:35
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
""" """
Module utils Module utils
@ -8,20 +29,159 @@ Module utils
""" """
import re
import numpy as np import numpy as np
from ase import Atoms, Atom from ase import Atom
from ase.visualize import view from ase import Atoms
class MsSpecAtoms(Atoms): from msspec.iodata import Data
def __init__(self, *args, **kwargs): from msspec.misc import LOGGER
Atoms.__init__(self, *args, **kwargs)
self.__absorber_index = None
def set_absorber(self, index):
self.__absorber_index = index
def get_absorber(self): class ForeignPotential(object):
return self.__absorber_index def __init__(self):
self.data = Data(title='Foreign Potential')
# Load exported potentials
# phagen_data is a dictionary with:
# self.phagen_data = {
# 'VintTotal' : <float>,
# 'VintCoulomb': <float>,
# 'RHOint' : <float>,
# 'types': [
# {
# 'Z' : <int>,
# 'RWS' : <float>,
# 'data': <np.array(..., 4, dtype=float)>
# },
# ...
# ...
# ...
# ]
# }
self.phagen_data = {'types': []}
def write(self, filename, prototypical_atoms):
LOGGER.debug(f"Writing Phagen input potential file: {filename}")
def append_atom_potential(atom):
Z = atom.number
# Find the right type (Z) in the phagen_data
itypes = []
for i, t in enumerate(self.phagen_data['types']):
if t['Z'] == Z:
itypes.append(i)
# Check now that we have only one type in the list
# otherwise we do not know yet how to deal with this.
assert len(itypes) > 0, f"Cannot find the data for atom with Z={Z}"
assert len(itypes) == 1, f"Too many datasets for atom with Z={Z}"
# So far so good, let's write the block
t = self.phagen_data['types'][itypes[0]]
s = "{:<7d}{:<10d}{:1.4f}\n".format(
t['Z'], len(t['data']), t['RWS'])
line_fmt = "{:+1.8e} " * 4 + "\n"
for row in t['data']:
s += line_fmt.format(*row)
return s
content = ""
# Start by writing the header line
content += "{:.2f} {:.4f} {:.2f}\n".format(
self.phagen_data['VintCoulomb'],
self.phagen_data['RHOint'],
self.phagen_data['VintTotal'])
# Then for each atom in the given prototypical cluster,
# write the data block
for atom in prototypical_atoms:
content += append_atom_potential(atom)
# Write the content to filename
try:
with open(filename, 'r') as fd:
old_content = fd.read()
except IOError:
old_content = ''
modified = False
if content != old_content:
with open(filename, 'w') as fd:
fd.write(content)
modified = True
return modified
class SPRKKRPotential(ForeignPotential):
def __init__(self, atoms, potfile, *exported_files):
super().__init__()
def read(content, pattern, types):
# compile the pattern for regex matching
pat = re.compile(pattern, re.MULTILINE)
# get the keys and data as list of strings
keys = re.split(r'\s+',
pat.search(content).group('KEYS').strip(' \n'))
txt = pat.search(content).group('DATA').strip('\n').split('\n')
data = []
for line in txt:
# Unpacking values
values_txt = re.split(r'\s+', line.strip())
data_dict = {}
for i, _ in enumerate(values_txt):
# Type casting
value = types[i].__call__(_)
data_dict[keys[i]] = value
# push to the list
data.append(data_dict)
return data
# load info in *.pot file
with open(potfile, 'r') as fd:
content = fd.read()
self.sites_data = read(content,
(r'^\s*SITES\s*\n((.*\n)+?\s*(?P<KEYS>IQ.*)\n'
r'(?P<DATA>(.*\n)+?))\*+'),
[int] + [float] * 3)
self.types_data = read(content,
(r'^\s*TYPES\s*\n(\s*(?P<KEYS>IT.*)\n'
r'(?P<DATA>(.*\n)+?))\*+'),
[int, str] + [int] * 4)
self.occ_data = read(content,
(r'^\s*OCCUPATION\s*\n(\s*(?P<KEYS>IQ.*)\n'
r'(?P<DATA>(.*\n)+?))\*+'),
[int] * 5 + [float])
for f in exported_files:
# get the IT from the filename
# m=re.match('SPRKKR-IT_(?P<IT>\d+)-PHAGEN.*', os.path.basename(f))
# it = int(m.group('IT'))
# load the content of the header (2 lines)
with open(f, 'r') as fd:
first_line = fd.readline().strip('\n')
second_line = fd.readline().strip('\n')
# load Coulomb and Total interstitial potential
pattern = (r'#\s*VMTZ_TOTAL\s*=\s*(?P<TOTAL>.*?)\s+'
r'VMTZ_Coulomb\s*=\s*(?P<COULOMB>.*?)\s+.*$')
m = re.match(pattern, first_line)
self.phagen_data.update(VintCoulomb=float(m.group('COULOMB')),
VintTotal=float(m.group('TOTAL')),
RHOint=0.)
# load Z, Wigner-Seitz radius from 2nd line of header
type_data = {}
_ = re.split(r'\s+', second_line.strip("# "))
type_data.update(Z=int(_[0]), RWS=float(_[3]))
# load the data
data = np.loadtxt(f, comments='#')
type_data.update(data=data)
self.phagen_data['types'].append(type_data)
class EmptySphere(Atom): class EmptySphere(Atom):
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
@ -42,7 +202,7 @@ def get_atom_index(atoms, x, y, z):
:rtype: int :rtype: int
""" """
# get all distances # get all distances
d = np.linalg.norm(atoms.get_positions() - np.array([x, y, z]), axis = 1) d = np.linalg.norm(atoms.get_positions() - np.array([x, y, z]), axis=1)
# get the index of the min distance # get the index of the min distance
i = np.argmin(d) i = np.argmin(d)
# return the index and the corresponding distance # return the index and the corresponding distance
@ -56,10 +216,9 @@ def center_cluster(atoms, invert=False):
with the origin being at the corner of the cell. with the origin being at the corner of the cell.
It is used in combination with the cut functions, which work only if the It is used in combination with the cut functions, which work only if the
origin is at the center of the cluster origin is at the center of the cluster
:param ase.Atoms atoms: an ASE Atoms object :param ase.Atoms atoms: an ASE Atoms object
:param bool invert: if True, performs the opposite translation (uncentering the cluster) :param bool invert: if True, performs the opposite translation
(uncentering the cluster)
""" """
for i, cell_vector in enumerate(atoms.get_cell()): for i, cell_vector in enumerate(atoms.get_cell()):
if invert: if invert:
@ -69,15 +228,6 @@ def center_cluster(atoms, invert=False):
def cut_sphere(atoms, radius, center=(0, 0, 0)): def cut_sphere(atoms, radius, center=(0, 0, 0)):
assert radius >= 0, "Please give a positive radius value"
radii = np.linalg.norm(atoms.positions - center, axis=1)
indices = np.where(radii <= radius)[0]
return atoms[indices]
def _cut_sphere(atoms, radius=None):
""" Removes all the atoms of an Atoms object outside a sphere with a """ Removes all the atoms of an Atoms object outside a sphere with a
given radius given radius
@ -87,28 +237,19 @@ def _cut_sphere(atoms, radius=None):
:return: The modified atom cluster :return: The modified atom cluster
:rtype: ase.Atoms :rtype: ase.Atoms
""" """
if radius is None: assert radius >= 0, "Please give a positive radius value"
raise ValueError("radius not set") radii = np.linalg.norm(atoms.positions - center, axis=1)
indices = np.where(radii <= radius)[0]
return atoms[indices]
new_atoms = atoms.copy()
del_list = []
for index, position in enumerate(new_atoms.positions):
if np.linalg.norm(position) > radius:
del_list.append(index)
del_list.reverse()
for index in del_list:
del new_atoms[index]
return new_atoms
def cut_cylinder(atoms, axis="z", radius=None): def cut_cylinder(atoms, axis="z", radius=None):
""" Removes all the atoms of an Atoms object outside a cylinder with a """ Removes all the atoms of an Atoms object outside a cylinder with a
given axis and radius given axis and radius
:param ase.Atoms atoms: an ASE Atoms object :param ase.Atoms atoms: an ASE Atoms object
:param str axis: string "x", "y", or "z". The axis of the cylinder, "z" by default :param str axis: string "x", "y", or "z". The axis of the cylinder,
"z" by default
:param float radius: the radius of the cylinder :param float radius: the radius of the cylinder
:return: The modified atom cluster :return: The modified atom cluster
@ -143,22 +284,23 @@ def cut_cylinder(atoms, axis="z", radius=None):
return new_atoms return new_atoms
def cut_cone(atoms, radius, z = 0):
def cut_cone(atoms, radius, z=0):
"""Shapes the cluster as a cone. """Shapes the cluster as a cone.
Keeps all the atoms of the input Atoms object inside a cone of based radius *radius* and of height *z*. Keeps all the atoms of the input Atoms object inside a cone of based
radius *radius* and of height *z*.
:param atoms: The cluster to modify. :param atoms: The cluster to modify.
:type atoms: :py:class:`ase.Atoms` :type atoms: :py:class:`ase.Atoms`
:param radius: The base cone radius in :math:`\mathring{A}`. :param radius: The base cone radius in :math:`\mathring{A}`. # noqa: W605
:type radius: float :type radius: float
:param z: The height of the cone in :math:`\mathring{A}`. :param z: The height of the cone in :math:`\mathring{A}`. # noqa: W605
:type z: float :type z: float
:return: A new cluster. :return: A new cluster.
:rtype: :py:class:`ase.Atoms` :rtype: :py:class:`ase.Atoms`
""" """
new_atoms = atoms.copy() new_atoms = atoms.copy()
origin = np.array((0, 0, 0))
max_theta = np.arctan(radius/(-z)) max_theta = np.arctan(radius/(-z))
u = np.array((0, 0, -z)) u = np.array((0, 0, -z))
@ -176,12 +318,12 @@ def cut_cone(atoms, radius, z = 0):
if theta <= max_theta: if theta <= max_theta:
indices.append(i) indices.append(i)
new_atoms = new_atoms[indices] new_atoms = new_atoms[indices]
new_atoms.translate(-u) # pylint: disable=invalid-unary-operand-type new_atoms.translate(-u)
return new_atoms return new_atoms
def cut_plane(atoms, x=None, y=None, z=None): def cut_plane(atoms, x=None, y=None, z=None):
""" Removes the atoms whose coordinates are higher (or lower, for a """ Removes the atoms whose coordinates are higher (or lower, for a
negative cutoff value) than the coordinates given for every dimension. negative cutoff value) than the coordinates given for every dimension.
@ -191,9 +333,9 @@ def cut_plane(atoms, x=None, y=None, z=None):
.. code-block:: python .. code-block:: python
cut_plane(atoms, x=[-5,5], y=3.6, z=0) cut_plane(atoms, x=[-5,5], y=3.6, z=0)
#every atom whose x-coordinate is higher than 5 or lower than -5, and/or # every atom whose x-coordinate is higher than 5 or lower than -5,
#y-coordinate is higher than 3.6, and/or z-coordinate is higher than 0 # and/or y-coordinate is higher than 3.6, and/or z-coordinate is higher
#is deleted. # than 0 is deleted.
:param ase.Atoms atoms: an ASE Atoms object :param ase.Atoms atoms: an ASE Atoms object
:param int x: x cutoff value :param int x: x cutoff value
@ -206,9 +348,11 @@ def cut_plane(atoms, x=None, y=None, z=None):
dim_names = ('x', 'y', 'z') dim_names = ('x', 'y', 'z')
dim_values = [x, y, z] dim_values = [x, y, z]
for i, (name, value) in enumerate(zip(dim_names, dim_values)): for i, (name, value) in enumerate(zip(dim_names, dim_values)):
assert isinstance(value, (int, float, list, tuple, type(None))), "Wrong type" assert isinstance(value, (int, float, list, tuple, type(None))), \
"Wrong type"
if isinstance(value, (tuple, list)): if isinstance(value, (tuple, list)):
assert len(value) == 2 and np.all([isinstance(el, (int, float, type(None))) for el in value]), \ assert len(value) == 2 and np.all(
[isinstance(el, (int, float, type(None))) for el in value]), \
"Wrong length" "Wrong length"
else: else:
try: try:
@ -216,7 +360,7 @@ def cut_plane(atoms, x=None, y=None, z=None):
dim_values[i] = [-np.inf, value] dim_values[i] = [-np.inf, value]
else: else:
dim_values[i] = [value, np.inf] dim_values[i] = [value, np.inf]
except: except Exception:
dim_values[i] = [value, value] dim_values[i] = [value, value]
if dim_values[i][0] is None: if dim_values[i][0] is None:
@ -227,15 +371,18 @@ def cut_plane(atoms, x=None, y=None, z=None):
dim_values = np.array(dim_values) dim_values = np.array(dim_values)
def constraint(coordinates): def constraint(coordinates):
return np.all(np.logical_and(coordinates >= dim_values[:,0], coordinates <= dim_values[:,1])) return np.all(np.logical_and(coordinates >= dim_values[:, 0],
coordinates <= dim_values[:, 1]))
indices = np.where(list(map(constraint, atoms.positions)))[0] indices = np.where(list(map(constraint, atoms.positions)))[0]
return atoms[indices] return atoms[indices]
def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
planes=0, shape='spherical'): planes=0, shape='spherical'):
"""Creates and returns a cluster based on an Atoms object and some parameters. """Creates and returns a cluster based on an Atoms object and some
parameters.
:param cluster: the Atoms object used to create the cluster :param cluster: the Atoms object used to create the cluster
:type cluster: Atoms object :type cluster: Atoms object
@ -245,7 +392,8 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
:type diameter: float :type diameter: float
:param planes: the number of planes of your cluster :param planes: the number of planes of your cluster
:type planes: integer :type planes: integer
:param emitter_plane: the plane where your emitter will be starting by 0 for the first plane :param emitter_plane: the plane where your emitter will be starting by 0
for the first plane
:type emitter_plane: integer :type emitter_plane: integer
See :ref:`hemispherical_cluster_faq` for more informations. See :ref:`hemispherical_cluster_faq` for more informations.
@ -255,7 +403,8 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
nmin = np.inf nmin = np.inf
for atom in cluster: for atom in cluster:
if ze - eps < atom.z < ze + eps and (atom.symbol == symbol or symbol == None): if ze - eps < atom.z < ze + eps and (atom.symbol == symbol or
symbol is None):
n = np.sqrt(atom.x**2 + atom.y**2) n = np.sqrt(atom.x**2 + atom.y**2)
if (n < nmin): if (n < nmin):
nmin = n nmin = n
@ -270,68 +419,115 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0,
eps = 0.01 # a useful small value eps = 0.01 # a useful small value
c = cell[:, 2].max() # a lattice parameter c = cell[:, 2].max() # a lattice parameter
a = cell[:, 0].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
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." # the number of planes in the cluster
p = np.alen(np.unique(np.round(cluster.get_positions()[:, 2], 4)))
# the symbol of your emitter
symbol = cluster[np.where(cluster.get_tags() == emitter_tag)[0][0]].symbol
assert (diameter != 0 or planes != 0), \
"At least one of diameter or planes parameter must be use."
if diameter == 0: if diameter == 0:
l = 1+2*(planes*c/p+1) # calculate the minimal diameter according to the number of planes # calculate the minimal diameter according to the number of planes
min_diameter = 1+2*(planes*c/p+1)
else: else:
l = diameter min_diameter = diameter
rep = int(2*l/min(a,c)) # number of repetition in each direction # number of repetition in each direction
cluster = cluster.repeat((rep, rep, rep)) # repeat the cluster rep = int(3*min_diameter/min(a, c))
center_cluster(cluster) # center the cluster # repeat the cluster
cluster.set_cell(cell) # reset the cell cluster = cluster.repeat((rep, rep, rep))
cluster = cut_plane(cluster, z=eps) # cut the cluster so that we have a centered surface
i = np.where(cluster.get_tags() == emitter_tag) # positions where atoms have the tag of the emitter_tag # center the cluster
all_ze = np.sort(np.unique(np.round(cluster.get_positions()[:, 2][i], 4))) # an array of all unique z corresponding to where we have the right atom's tag center_cluster(cluster)
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # an array of all unique z
n = np.where(all_z == all_z.max())[0][0] - np.where(all_z == all_ze.max())[0][0] # calculate the number of planes above the emitter's plane # reset the cell
ze = all_ze.max() # the height of the emitter's plane cluster.set_cell(cell)
# cut the cluster so that we have a centered surface
cluster = cut_plane(cluster, z=eps)
# if the number of planes above the emitter's plane is smaller than it must be, recalculate n and ze # positions where atoms have the tag of the emitter_tag
while n < emitter_plane: i = np.where(cluster.get_tags() == emitter_tag)
all_ze = all_ze[:-1]
n = np.where(all_z == all_z.max())[0][0] - np.where(all_z == all_ze.max())[0][0] # an array of all unique z corresponding to where we have the right
# atom's tag
all_ze = np.sort(np.unique(np.round(cluster.get_positions()[:, 2][i], 4)))
# an array of all unique z
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))
# calculate the number of planes above the emitter's plane
n = np.where(all_z == all_z.max())[0][0] - np.where(
all_z == all_ze.max())[0][0]
# the height of the emitter's plane
ze = all_ze.max() ze = all_ze.max()
# if the number of planes above the emitter's plane is smaller than it must
# be, recalculate n and ze
while n < emitter_plane:
all_ze = all_ze[:-1]
n = np.where(all_z == all_z.max())[0][0] - np.where(
all_z == all_ze.max())[0][0]
ze = all_ze.max()
tx, ty = get_xypos(cluster, ze, symbol) # values of x and y of the emitter # values of x and y of the emitter
Atoms.translate(cluster, [-tx, -ty, 0]) # center the cluster on the emitter tx, ty = get_xypos(cluster, ze, symbol)
z_cut = all_z[np.where(all_z == all_ze.max())[0][0] + emitter_plane] # calculate where to cut to get the right number of planes above the emitter # center the cluster on the emitter
Atoms.translate(cluster, [0, 0, -z_cut]) # translate the surface at z=0 Atoms.translate(cluster, [-tx, -ty, 0])
cluster = cut_plane(cluster, z=eps) # cut the planes above those we want to keep
# calculate where to cut to get the right number of planes above the
# emitter
z_cut = all_z[np.where(all_z == all_ze.max())[0][0] + emitter_plane]
# translate the surface at z=0
Atoms.translate(cluster, [0, 0, -z_cut])
# cut the planes above those we want to keep
cluster = cut_plane(cluster, z=eps)
radius = diameter/2 radius = diameter/2
if planes!=0: if planes != 0:
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # an array of all unique remaining z # an array of all unique remaining z
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))
zplan = all_z[-planes] zplan = all_z[-planes]
xplan, yplan = get_xypos(cluster, zplan) xplan, yplan = get_xypos(cluster, zplan)
radius = np.sqrt(xplan**2 + yplan**2 + zplan**2) radius = np.sqrt(xplan**2 + yplan**2 + zplan**2)
if diameter!=0:
assert (radius <= diameter/2), "The number of planes is too high compared to the diameter." if diameter != 0:
assert (radius <= diameter/2), ("The number of planes is too high "
"compared to the diameter.")
radius = max(radius, diameter/2) radius = max(radius, diameter/2)
if shape in ('spherical'): if shape in ('spherical'):
cluster = cut_sphere(cluster, radius=radius + eps) # cut a sphere in our cluster with the diameter which is indicate in the parameters # cut a sphere in our cluster with the diameter which is indicate in
# the parameters
cluster = cut_sphere(cluster, radius=radius + eps)
elif shape in ('cylindrical'): elif shape in ('cylindrical'):
cluster = cut_cylinder(cluster, radius=radius + eps) # cut a sphere in our cluster with the diameter which is indicate in the parameters # cut a sphere in our cluster with the diameter which is indicate in
# the parameters
cluster = cut_cylinder(cluster, radius=radius + eps)
else: else:
raise NameError('Unkknown shape specifier: \"{}\"'.format(shape)) raise NameError('Unkknown shape specifier: \"{}\"'.format(shape))
if planes!=0: if planes != 0:
zcut = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))[::-1][planes-1] - eps # calculate where to cut to get the right number of planes # calculate where to cut to get the right number of planes
cluster = cut_plane(cluster, z=zcut) # cut the right number of planes positions = np.unique(np.round(cluster.get_positions()[:, 2], 4))
zcut = np.sort(positions)[::-1][planes-1] - eps
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # an array of all unique remaining z # cut the right number of planes
assert emitter_plane < np.alen(all_z), "There are not enough existing plans." cluster = cut_plane(cluster, z=zcut)
# an array of all unique remaining z
all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4)))
assert emitter_plane < np.alen(all_z), ("There are not enough existing "
"plans.")
ze = all_z[- emitter_plane - 1] # the z-coordinate of the emitter 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) Atoms.translate(cluster, [0, 0, -ze]) # put the emitter in (0,0,0)

31
src/msspec/version.py
View File

@ -1,10 +1,33 @@
# coding: utf-8 #!/usr/bin/env python
# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a: #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# msspec is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this msspec. If not, see <http://www.gnu.org/licenses/>.
#
# Source file : src/msspec/version.py
# Last modified: ven. 10 avril 2020 17:34:38
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
from setuptools_scm import get_version
from pkg_resources import parse_version, DistributionNotFound, get_distribution
import os import os
from pkg_resources import DistributionNotFound
from pkg_resources import get_distribution
from pkg_resources import parse_version
from setuptools_scm import get_version
# find the version number # find the version number
# 1- Try to read it from the git info # 1- Try to read it from the git info
# 2- If it fails, try to read it from the distribution file # 2- If it fails, try to read it from the distribution file

39
src/setup.py
View File

@ -1,22 +1,31 @@
# coding: utf-8 #!/usr/bin/env python
# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : #
# Copyright © 2016-2020 - Rennes Physics Institute
#
# This file is part of msspec.
#
# msspec is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
from setuptools import setup, Extension, find_packages # msspec is distributed in the hope that it will be useful,
from distutils.file_util import copy_file # but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
from distutils.command.build import build as _build # You should have received a copy of the GNU General Public License
from setuptools.command.build_ext import build_ext as _build_ext # along with this msspec. If not, see <http://www.gnu.org/licenses/>.
from setuptools.command.install import install as _install #
from distutils.command.clean import clean as _clean # Source file : src/setup.py
# Last modified: mar. 07 avril 2020 17:01:42
# Committed by : "Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>"
import subprocess from setuptools import setup, find_packages
import traceback from version import __version__
import os
import sys import sys
import glob
sys.path.insert(0, "msspec") sys.path.insert(0, "msspec")
from version import __version__
with open('setup_requirements.txt', 'r') as fd: with open('setup_requirements.txt', 'r') as fd:
SETUP_REQUIREMENTS = fd.read().strip().split('\n') SETUP_REQUIREMENTS = fd.read().strip().split('\n')
@ -37,8 +46,8 @@ if __name__ == "__main__":
maintainer='Sylvain Tricot', maintainer='Sylvain Tricot',
maintainer_email='sylvain.tricot@univ-rennes1.fr', maintainer_email='sylvain.tricot@univ-rennes1.fr',
url='https://msspec.cnrs.fr', url='https://msspec.cnrs.fr',
description=('A multiple scattering package for sepectroscopies using ' description=('A multiple scattering package for sepectroscopies '
'electrons to probe materials'), 'using electrons to probe materials'),
long_description="""MsSpec is a Fortran package to compute the long_description="""MsSpec is a Fortran package to compute the
cross-section of several spectroscopies involving one (or more) cross-section of several spectroscopies involving one (or more)
electron(s) as the probe. This package provides a python interface to electron(s) as the probe. This package provides a python interface to