diff --git a/doc/source/tutorials/copper/log.log b/doc/source/tutorials/copper/log.log deleted file mode 100644 index 9f405a2..0000000 --- a/doc/source/tutorials/copper/log.log +++ /dev/null @@ -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 ** - - diff --git a/src/msspec/__init__.py b/src/msspec/__init__.py index fe5175a..9a852c4 100644 --- a/src/msspec/__init__.py +++ b/src/msspec/__init__.py @@ -1,8 +1,30 @@ -import ase -from .version import __version__ -from .misc import LOGGER +#!/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 . +# +# Source file : src/msspec/__init__.py +# Last modified: ven. 10 avril 2020 17:22:12 +# Committed by : "Sylvain Tricot " + + +import ase + +from msspec.misc import LOGGER +from msspec.version import __version__ -__sha__ = "$Id$" def init_msspec(): LOGGER.debug('Initialization of the msspec module') diff --git a/src/msspec/calcio.py b/src/msspec/calcio.py index ce14e09..324712a 100644 --- a/src/msspec/calcio.py +++ b/src/msspec/calcio.py @@ -1,4 +1,5 @@ # coding: utf-8 +# vim: set et ts=4 sw=4 fdm=indent mouse=a cc=+1 tw=80: """ Module calcio @@ -6,13 +7,16 @@ Module calcio """ -import numpy as np from datetime import datetime -import ase 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 + class PhagenIO(object): def __init__(self, phagen_parameters, malloc_parameters): self.parameters = phagen_parameters @@ -27,8 +31,9 @@ class PhagenIO(object): self.nlmax = 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 + LOGGER.debug("Writing Phagen input file") atoms = self.parameters.atoms if not atoms.has('mt_radii'): for a in atoms: @@ -58,7 +63,7 @@ class PhagenIO(object): # other sections continue 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)) else: s = ' {}={:s}'.format(name, str(parameter)) @@ -167,7 +172,8 @@ class PhagenIO(object): lines = pattern.split(content) # 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.ne = ne self.ipot = ipot @@ -181,7 +187,7 @@ class PhagenIO(object): if not re.match(r'^\d+$', n): numbers.append(float(n)) 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) # construct the data array @@ -265,8 +271,27 @@ class PhagenIO(object): proto_index = int(data[i, -1]) proto_indices.append(proto_index) 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( - 1, self.nat + 1)])) + self.nateqm = int(np.max([ + 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'): a_index = 0 @@ -279,11 +304,13 @@ class PhagenIO(object): a_index += 1 d = d.split() 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': []} pot_data.append(a) 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 for _pot_data in pot_data: @@ -292,27 +319,30 @@ class PhagenIO(object): return pot_data + class SpecIO(object): def __init__(self, parameters, malloc_parameters, phagenio): self.parameters = parameters self.malloc_parameters = malloc_parameters 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): if center: - s = ('{}*{:^%ds}*\n' % (width - shift - 2)).format(' ' * shift, t) + s = ('{}*{:^%ds}*\n' % (width - shift - 2) + ).format(' ' * shift, t) else: - s = ('{}*{:%ds}*\n' % (width - shift - 2)).format(' ' * shift, t) + s = ('{}*{:%ds}*\n' % (width - shift - 2) + ).format(' ' * shift, t) return s def rule(tabs=(5, 10, 10, 10, 10), symbol='=', shift=4, width=79): s = ' ' * shift + '*' + symbol * (width - shift - 2) + '*\n' t = np.cumsum(tabs) + shift - l = list(s) + sep = list(s) for i in t: - l[i] = '+' - return ''.join(l) + sep[i] = '+' + return ''.join(sep) def fillstr(a, b, index, justify='left'): alist = list(a) @@ -426,7 +456,8 @@ class SpecIO(object): content += rule() line = create_line("SPECTRO,ISPIN,IDICHR,IPOL") 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_ipol), 39, 'left') 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_imoy')), 19, '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 += rule() - content += title(' ' * 22 + 'LEED EXPERIMENTAL PARAMETERS :', center=False) + content += title(' ' * 22 + 'LEED EXPERIMENTAL PARAMETERS :', + center=False) content += rule() line = create_line("IPHI,ITHETA,IE,IFTHET") line = fillstr(line, str(p.get_parameter('leed_iphi')), 9, 'left') @@ -496,48 +529,57 @@ class SpecIO(object): content += line 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_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_r0')), 39, 'decimal') content += line 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_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_r1')), 39, 'decimal') content += line 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_phiini')), 19, 'decimal') + line = fillstr(line, + str(p.get_parameter('leed_phiini')), 19, 'decimal') content += line 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_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, 'decimal') content += line content += rule() - content += title(' ' * 21 + 'EXAFS EXPERIMENTAL PARAMETERS :', center=False) + content += title(' ' * 21 + 'EXAFS EXPERIMENTAL PARAMETERS :', + center=False) content += rule() 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_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, 'decimal') content += line 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_ekini')), 19, 'decimal') - line = fillstr(line, str(p.get_parameter('exafs_ekfin')), 29, 'decimal') + line = fillstr(line, + 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, 'decimal') content += line content += rule() - content += title(' ' * 22 + 'AED EXPERIMENTAL PARAMETERS :', center=False) + content += title(' ' * 22 + 'AED EXPERIMENTAL PARAMETERS :', + center=False) content += rule() line = create_line("EDGE_C,EDGE_I,EDGE_A") 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_method')), 19, 'left') 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 line = create_line("N_MAX,N_ITER,N_TABLE,SHIFT") line = fillstr(line, str(p.get_parameter('eigval_nmax')), 9, 'left') @@ -673,7 +716,8 @@ class SpecIO(object): thfwd_arr = np.ones((nat)) path_filtering = p.extra_parameters['calculation'].get_parameter( '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) else: ibwd_arr = np.zeros((nat), dtype=np.int) @@ -694,7 +738,8 @@ class SpecIO(object): 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_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') content += line line = create_line("ILENGTH,RLENGTH,UNLENGTH") @@ -713,7 +758,8 @@ class SpecIO(object): idcm = format(2, 'd') temp = format(0., '.2f') td = format(500., '.2f') - LOGGER.warning('Vibrational damping is disabled for this calculation.') + LOGGER.warning('Vibrational damping is disabled for this ' + 'calculation.') else: idwsph = str(p.get_parameter('calc_idwsph')) idcm = str(p.get_parameter('calc_idcm')) @@ -757,8 +803,9 @@ class SpecIO(object): content += title(' ' * 17 + 'INPUT FILES (PHD, EXAFS, LEED, AED, ' 'APECS) :', center=False) content += rule(tabs=(), symbol='-') - content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', - center=False) + content += title( + ' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', + center=False) content += rule(tabs=(5, 23, 7, 10)) line = create_line("DATA FILE,UNIT") line = fillstr(line, str(p.get_parameter('input_data')), 9, 'right') @@ -791,8 +838,9 @@ class SpecIO(object): center=False) content += title(' ' * 28 + '(AUGER ELECTRON)', center=False) content += rule(tabs=(), symbol='-') - content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', - center=False) + content += title( + ' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', + center=False) content += rule(tabs=(5, 23, 7, 10)) line = create_line("PHASE SHIFTS/TL FILE,UNIT") line = fillstr(line, str(p.get_parameter('input2_tl')), 9, 'right') @@ -810,8 +858,9 @@ class SpecIO(object): content += title(' ' * 29 + 'OUTPUT FILES :', center=False) content += rule(tabs=(), symbol='-') - content += title(' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', - center=False) + content += title( + ' ' * 8 + 'NAME' + ' ' * 20 + 'UNIT' + ' ' * 16 + 'TYPE', + center=False) content += rule(tabs=(5, 23, 7, 10)) line = create_line("CONTROL FILE,UNIT") line = fillstr(line, str(p.get_parameter('output_log')), 9, 'right') @@ -826,7 +875,8 @@ class SpecIO(object): line = fillstr(line, str(p.get_parameter('output_unit11')), 39, 'left') content += line 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') content += line content += rule(tabs=(5, 23, 7, 10)) @@ -874,15 +924,13 @@ class SpecIO(object): # backup the content in memory old_content = content - """ - for key in ('NATP_M', 'NATCLU_M', 'NE_M', 'NEMET_M', 'LI_M', 'NL_M', - 'NO_ST_M'): - required = requirements[key] - limit = self.malloc_parameters.get_parameter(key).value - value = required if required > limit else limit - content = re.sub(r'({:s}\s*=\s*)\d+'.format(key), - r'\g<1>{:d}'.format(value), content) - """ + # for key in ('NATP_M', 'NATCLU_M', 'NE_M', 'NEMET_M', 'LI_M', 'NL_M', + # 'NO_ST_M'): + # required = requirements[key] + # limit = self.malloc_parameters.get_parameter(key).value + # value = required if required > limit else limit + # content = re.sub(r'({:s}\s*=\s*)\d+'.format(key), + # r'\g<1>{:d}'.format(value), content) 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', @@ -895,7 +943,6 @@ class SpecIO(object): content = re.sub(r'({:s}\s*=\s*)\d+'.format(key), r'\g<1>{:d}'.format(value), content) - modified = False if content != old_content: with open(filename, 'w') as fd: @@ -945,7 +992,7 @@ class SpecIO(object): data = np.loadtxt(filename, skiprows=skip, unpack=True) if len(data.shape) <= 1: - data = data.reshape((2,-1)) + data = data.reshape((2, -1)) return data 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+)') with open(filename, 'r') as fd: content = fd.read() - #return pat.findall(content.replace('\n', '__cr__')) return pat.findall(content) - - - diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 8d00198..e0c2188 100644 --- a/src/msspec/calculator.py +++ b/src/msspec/calculator.py @@ -1,5 +1,25 @@ -# coding: utf-8 -# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a tw=80: +#!/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 . +# +# Source file : src/msspec/calculator.py +# Last modified: ven. 10 avril 2020 17:19:24 +# Committed by : "Sylvain Tricot " + + """ Module calculator ================= @@ -21,54 +41,58 @@ For more information on MsSpec, follow this """ - - +import inspect import os +import re import shutil import sys -import re -import inspect - -from subprocess import Popen, PIPE -from shutil import copyfile, rmtree -from datetime import datetime import time -from io import StringIO 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.atoms - +import ase.data import numpy as np - +from ase.calculators.calculator import Calculator +from terminaltables.ascii_table import AsciiTable 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.misc import (UREG, LOGGER, get_call_info, - get_level_from_electron_configuration, - XRaySource, set_log_output, log_process_output) -from msspec.utils import get_atom_index - -from msspec.parameters import (PhagenParameters, PhagenMallocParameters, - SpecParameters, SpecMallocParameters, - GlobalParameters, MuffintinParameters, - TMatrixParameters, SourceParameters, - DetectorParameters, ScanParameters, - CalculationParameters, - PEDParameters, EIGParameters) -from msspec.calcio import PhagenIO, SpecIO - +from msspec.data import electron_be +from msspec.misc import get_call_info +from msspec.misc import get_level_from_electron_configuration +from msspec.misc import log_process_output +from msspec.misc import LOGGER +from msspec.misc import set_log_output +from msspec.misc import UREG +from msspec.misc import XRaySource +from msspec.parameters import CalculationParameters +from msspec.parameters import DetectorParameters +from msspec.parameters import EIGParameters +from msspec.parameters import GlobalParameters +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.spec.fortran import _phd_se_noso_nosp_nosym, _phd_mi_noso_nosp_nosym from msspec.spec.fortran import _eig_mi from msspec.spec.fortran import _eig_pw - - -from terminaltables.ascii_table import AsciiTable +from msspec.spec.fortran import _phd_mi_noso_nosp_nosym +from msspec.spec.fortran import _phd_se_noso_nosp_nosym +from msspec.utils import get_atom_index class _MSCALCULATOR(Calculator): @@ -78,9 +102,9 @@ class _MSCALCULATOR(Calculator): """ implemented_properties = ['', ] __data = {} - + def __init__(self, spectroscopy='PED', algorithm='expansion', - polarization=None, dichroism=None, spinpol=False, + polarization=None, dichroism=None, spinpol=False, folder='./calc', txt='-', **kwargs): stdout = sys.stdout if isinstance(txt, str) and txt != '-': @@ -96,7 +120,7 @@ class _MSCALCULATOR(Calculator): ######################################################################## # Init the upper class Calculator.__init__(self, **kwargs) - + ######################################################################## LOGGER.debug(' create low level parameters') ######################################################################## @@ -104,15 +128,15 @@ class _MSCALCULATOR(Calculator): self.phagen_malloc_parameters = PhagenMallocParameters() self.spec_parameters = SpecParameters() self.spec_malloc_parameters = SpecMallocParameters() - + ######################################################################## LOGGER.debug(' create higher level parameters') ######################################################################## self.tmatrix_parameters = TMatrixParameters(self.phagen_parameters) self.muffintin_parameters = MuffintinParameters(self.phagen_parameters, self.spec_parameters) - - + + self.global_parameters = GlobalParameters(self.phagen_parameters, self.spec_parameters) @@ -124,25 +148,25 @@ class _MSCALCULATOR(Calculator): self.spec_parameters) else: raise NameError('No such spectroscopy') - + self.source_parameters = SourceParameters(self.global_parameters, self.phagen_parameters, self.spec_parameters) - + self.detector_parameters = DetectorParameters(self.global_parameters, self.phagen_parameters, self.spec_parameters) - + self.scan_parameters = ScanParameters(self.global_parameters, self.phagen_parameters, - self.spec_parameters) - + self.spec_parameters) + self.calculation_parameters = CalculationParameters( self.global_parameters, self.phagen_parameters, self.spec_parameters) # initialize all parameters with defaults values LOGGER.info("Set default values =========================================") - for p in (list(self.global_parameters) + + for p in (list(self.global_parameters) + list(self.muffintin_parameters) + list(self.tmatrix_parameters) + list(self.source_parameters) + @@ -160,15 +184,15 @@ class _MSCALCULATOR(Calculator): self.global_parameters.dichroism = dichroism self.global_parameters.spinpol = spinpol self.global_parameters.folder = folder - - - - self.phagenio = PhagenIO(self.phagen_parameters, + + + + self.phagenio = PhagenIO(self.phagen_parameters, self.phagen_malloc_parameters) - self.specio = SpecIO(self.spec_parameters, + self.specio = SpecIO(self.spec_parameters, self.spec_malloc_parameters, self.phagenio) - + ######################################################################## LOGGER.debug(' create a space dedicated to the calculation') ######################################################################## @@ -180,7 +204,7 @@ class _MSCALCULATOR(Calculator): os.makedirs(self.tmp_folder) os.makedirs(os.path.join(self.tmp_folder, 'input')) os.makedirs(os.path.join(self.tmp_folder, 'output')) - #copyfile(os.path.join(self.msspec_folder, 'ase', 'Makefile'), + #copyfile(os.path.join(self.msspec_folder, 'ase', 'Makefile'), # os.path.join(self.tmp_folder, 'Makefile')) #os.chdir(self.tmp_folder) @@ -202,7 +226,7 @@ class _MSCALCULATOR(Calculator): dichro_spinpol = True #spin = 'NO' - #if spinpol or dichro_spinpol: + #if spinpol or dichro_spinpol: # spin = 'YES' #if spin == 'YES': @@ -225,7 +249,7 @@ class _MSCALCULATOR(Calculator): ######################################################################## LOGGER.debug(' initialization done.\n') ######################################################################## - + def _make(self, target): LOGGER.debug(get_call_info(inspect.currentframe())) os.chdir(self.tmp_folder) @@ -252,11 +276,11 @@ class _MSCALCULATOR(Calculator): LOGGER.error("Unable to successfully run the target: {}".format(target)) sys.exit(1) - + def _guess_ke(self, level): """ Try to guess the kinetic energy based on the level and the source energy. If the kinetic energy cannot be infered - because the level is not reported in the database, the returned + because the level is not reported in the database, the returned value is None. """ try: @@ -274,39 +298,38 @@ class _MSCALCULATOR(Calculator): ke = source_energy - binding_energy - wf #return np.array(ke, dtype=np.float).flatten() return ke - - + + 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') - - #mod0 = self.phagenio.write_include_file(filename=include_fname) - mod1 = self.phagenio.write_input_file(filename=input_fname) - - self.modified = self.modified or mod1 #or mod0 or mod1 - - #if self.modified: - if mod1: - # run phagen - #self._make('tmatrix') + modified = self.phagenio.write_input_file(filename=input_fname) + + # Write the external potential file + 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) + + # Run if input files changed + if modified: + # Run phagen os.chdir(os.path.join(self.tmp_folder, 'output')) - # copy external potential file to the workspace - if self.tmatrix_parameters.potential in ('lmto', 'spkkr', 'msf'): - fname = os.path.join(self.init_folder, - self.tmatrix_parameters.potential_file) - shutil.copy(fname, 'fort.2') do_phagen() - # rename some output files to be more explicit - #os.rename('fort.10', 'cluster.clu') - #os.rename('fort.35', 'tmatrix.tl') - #os.rename('fort.55', 'tmatrix.rad') - #shutil.copy('fort.10', 'cluster.clu') + # Copy some output files to be more explicit shutil.copy('clus/clus.out', 'cluster.clu') shutil.copy('fort.35', 'tmatrix.tl') 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) - + def run_spec(self, malloc={}): def get_li(level): @@ -317,14 +340,14 @@ class _MSCALCULATOR(Calculator): #include_fname = os.path.join(self.tmp_folder, 'src/spec.inc') input_fname = os.path.join(self.tmp_folder, 'input/spec.dat') kdirs_fname = os.path.join(self.tmp_folder, 'input/kdirs.dat') - + mod0 = self.specio.write_input_file(filename=input_fname) #mod1 = self.specio.write_include_file(filename=include_fname) mod2 = self.specio.write_kdirs_file(filename=kdirs_fname) - + #self.modified = self.modified or mod0 or mod1 or mod2 self.modified = self.modified or mod0 or mod2 - + #self._make('tmatrix') #self._make('bin/spec') #t0 = time.time() @@ -424,25 +447,29 @@ class _MSCALCULATOR(Calculator): def get_tmatrix(self): LOGGER.info("Getting the TMatrix...") 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() - - filename = os.path.join(self.tmp_folder, 'output/tmatrix.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 = self.phagenio.tl tl_threshold = self.tmatrix_parameters.get_parameter('tl_threshold') if tl_threshold.value != None: - LOGGER.debug(" applying tl_threshold to %s...", + LOGGER.debug(" applying tl_threshold to %s...", tl_threshold.value) - go_on = True + go_on = True while go_on: go_on = False for ia in range(self.phagenio.nat): - for ie in range(self.phagenio.ne): + for ie in range(self.phagenio.ne): last_tl = tl[ia][ie][-1, -2:] # convert to complex last_tl = last_tl[0] + last_tl[1] * 1j @@ -450,36 +477,36 @@ class _MSCALCULATOR(Calculator): # remove last line of tl tl[ia][ie] = tl[ia][ie][:-1, :] go_on = True - + max_tl = self.tmatrix_parameters.get_parameter('max_tl').value cluster = self.phagen_parameters.get_parameter('atoms').value proto_indices = cluster.get_array('proto_indices') - + if max_tl != None: LOGGER.debug(" applying max_tl: %s", max_tl) for ia in range(self.phagenio.nat): for ie in range(self.phagenio.ne): try: - # for each set of tl: + # for each set of tl: # 1. get the symbol of the prototipical atom j = np.where(proto_indices == ia+1)[0] - symbol = cluster[j][0].symbol + symbol = cluster[j][0].symbol # 2. get the number of max tl allowed ntl = max_tl[symbol] # 3. reshape the tl set accordingly tl[ia][ie] = tl[ia][ie][:ntl, :] except KeyError: - pass + pass self.phagenio.write_tl_file( os.path.join(self.tmp_folder, 'output/tmatrix.tl')) - + # update spec extra parameters here self.spec_parameters.set_parameter('extra_nat', self.phagenio.nat) self.spec_parameters.set_parameter('extra_nlmax', self.phagenio.nlmax) - - + + def set_atoms(self, atoms): """Defines the cluster on which the calculator will work. @@ -502,7 +529,7 @@ class _MSCALCULATOR(Calculator): """ _ = [] - for section in ('global', 'muffintin', 'tmatrix', 'spectroscopy', + for section in ('global', 'muffintin', 'tmatrix', 'spectroscopy', 'source', 'detector', 'scan', 'calculation'): parameters = getattr(self, section + '_parameters') for p in parameters: @@ -514,7 +541,7 @@ class _MSCALCULATOR(Calculator): self.atoms.info['absorber'] = self.atoms.absorber self.atoms.write(clusbuf, format='xyz') dset.add_parameter(group='Cluster', name='cluster', value=clusbuf.getvalue(), hidden="True") - + def shutdown(self): """Removes the temporary folder and all its content. @@ -539,7 +566,7 @@ class _PED(_MSCALCULATOR): """This class creates a calculator object for PhotoElectron DIffraction spectroscopy. - :param algorithm: The algorithm to use for the computation. See + :param algorithm: The algorithm to use for the computation. See :ref:`globalparameters-algorithm` for more details about the allowed values and the type. @@ -573,24 +600,24 @@ class _PED(_MSCALCULATOR): spinpol=spinpol, folder=folder, txt=txt) self.iodata = iodata.Data('PED Simulation') - - def _get_scan(self, scan_type='theta', phi=0, - theta=np.linspace(-70, 70, 141), level=None, + + def _get_scan(self, scan_type='theta', phi=0, + theta=np.linspace(-70, 70, 141), level=None, kinetic_energy=None, data=None, malloc={}): LOGGER.info("Computting the %s scan...", scan_type) if data: self.iodata = data - + if kinetic_energy is None: # try to guess the kinetic energy kinetic_energy = self._guess_ke(level) - + # if still None... if kinetic_energy is None: LOGGER.error('Unable to guess the kinetic energy!') raise ValueError('You must define a kinetic_energy value.') - + # update the parameters self.scan_parameters.set_parameter('kinetic_energy', kinetic_energy) all_ke = self.scan_parameters.get_parameter('ke_array') @@ -599,29 +626,29 @@ class _PED(_MSCALCULATOR): raise ValueError('Kinetic energy is < 0! ({})'.format( kinetic_energy)) self.scan_parameters.set_parameter('type', scan_type) - + # make sure there is only one energy point in scatf scan if scan_type == 'scatf': assert len(all_ke) == 1, ('kinetic_energy should not be an array ' 'in scatf scan') - - + + if scan_type != 'scatf': self.scan_parameters.set_parameter('phi', phi) self.scan_parameters.set_parameter('theta', theta) - + self.spectroscopy_parameters.set_parameter('level', level) - + self.get_tmatrix() self.run_spec(malloc) - + # Now load the data ndset = len(self.iodata) dset = self.iodata.add_dset('{} scan [{:d}]'.format(scan_type, ndset)) for p in self.get_parameters(): - bundle = {'group': str(p.group), + bundle = {'group': str(p.group), 'name': str(p.name), - 'value': str(p.value), + 'value': str(p.value), 'unit': '' if p.unit is None else str(p.unit)} dset.add_parameter(**bundle) if scan_type in ('theta', 'phi', 'energy'): @@ -636,7 +663,7 @@ class _PED(_MSCALCULATOR): data = data[:, [1, 4, 5, 6, 8]].T _proto, _sf_real, _sf_imag, _theta, _energy = data _sf = _sf_real + _sf_imag * 1j - dset.add_columns(proto_index=_proto, sf_real=np.real(_sf), + dset.add_columns(proto_index=_proto, sf_real=np.real(_sf), sf_imag=np.imag(_sf), sf_module=np.abs(_sf), theta=_theta, energy=_energy) elif scan_type in ('theta_phi',): @@ -645,7 +672,7 @@ class _PED(_MSCALCULATOR): #theta_c, phi_c = data[[2, 3], :] #xsec_c = data[-1, :] #dirsig_c = data[-2, :] - + #dset.add_columns(theta=theta_c) #dset.add_columns(phi=phi_c) #dset.add_columns(cross_section=xsec_c) @@ -665,14 +692,14 @@ class _PED(_MSCALCULATOR): xlabel = r'Angle $\theta$($\degree$)' ylabel = r'Signal (a. u.)' - view = dset.add_view("E = {:.2f} eV".format(ke), title=title, - xlabel=xlabel, ylabel=ylabel) + view = dset.add_view("E = {:.2f} eV".format(ke), title=title, + xlabel=xlabel, ylabel=ylabel) for angle_phi in self.scan_parameters.get_parameter( 'phi').value: where = ("energy=={:.2f} and phi=={:.2f}" "").format(ke, angle_phi) legend = r'$\phi$ = {:.1f} $\degree$'.format(angle_phi) - view.select('theta', 'cross_section', where=where, + view.select('theta', 'cross_section', where=where, legend=legend) if scan_type == 'phi': absorber_symbol = self.atoms[self.atoms.absorber].symbol @@ -681,16 +708,16 @@ class _PED(_MSCALCULATOR): xlabel = r'Angle $\phi$($\degree$)' ylabel = r'Signal (a. u.)' - view = dset.add_view("E = {:.2f} eV".format(ke), title=title, - xlabel=xlabel, ylabel=ylabel) + view = dset.add_view("E = {:.2f} eV".format(ke), title=title, + xlabel=xlabel, ylabel=ylabel) for angle_theta in self.scan_parameters.get_parameter( 'theta').value: where = ("energy=={:.2f} and theta=={:.2f}" "").format(ke, angle_theta) legend = r'$\theta$ = {:.1f} $\degree$'.format(angle_theta) - view.select('phi', 'cross_section', where=where, + view.select('phi', 'cross_section', where=where, legend=legend) - + if scan_type == 'theta_phi': absorber_symbol = self.atoms[self.atoms.absorber].symbol title = ('Stereographic projection of {}({}) at {:.2f} eV' @@ -698,20 +725,20 @@ class _PED(_MSCALCULATOR): xlabel = r'Angle $\phi$($\degree$)' ylabel = r'Signal (a. u.)' - view = dset.add_view("E = {:.2f} eV".format(ke), title=title, - xlabel=xlabel, ylabel=ylabel, + view = dset.add_view("E = {:.2f} eV".format(ke), title=title, + xlabel=xlabel, ylabel=ylabel, projection='stereo', colorbar=True, autoscale=True) view.select('theta', 'phi', 'cross_section') - - + + if scan_type == 'scatf': for i in range(self.phagenio.nat): proto_index = i+1 title = 'Scattering factor at {:.3f} eV'.format(kinetic_energy) - view = dset.add_view("Proto. atom #{:d}".format(proto_index), + view = dset.add_view("Proto. atom #{:d}".format(proto_index), title=title, projection='polar') - where = "proto_index=={:d}".format(proto_index) + where = "proto_index=={:d}".format(proto_index) view.select('theta', 'sf_module', where=where, legend=r'$|f(\theta)|$') view.select('theta', 'sf_real', where=where, @@ -795,13 +822,13 @@ class _PED(_MSCALCULATOR): return self.iodata - def get_scattering_factors(self, level='1s', kinetic_energy=None, + def get_scattering_factors(self, level='1s', kinetic_energy=None, data=None): """Computes the scattering factors of all prototypical atoms in the cluster. This function computes the real and imaginery parts of the scattering - factor as well as its modulus for each non symetrically equivalent atom + factor as well as its modulus for each non symetrically equivalent atom in the cluster. The results are stored in the *data* object if provided as a parameter. @@ -814,11 +841,11 @@ class _PED(_MSCALCULATOR): argument or a new :py:class:`iodata.Data` object. """ - data = self._get_scan(scan_type='scatf', level=level, data=data, + data = self._get_scan(scan_type='scatf', level=level, data=data, kinetic_energy=kinetic_energy) return data - - def get_theta_scan(self, phi=0, theta=np.linspace(-70, 70, 141), + + def get_theta_scan(self, phi=0, theta=np.linspace(-70, 70, 141), level=None, kinetic_energy=None, data=None): """Computes a polar scan of the emitted photoelectrons. @@ -838,7 +865,7 @@ class _PED(_MSCALCULATOR): data = self._get_scan(scan_type='theta', level=level, theta=theta, phi=phi, kinetic_energy=kinetic_energy, data=data) return data - + def get_phi_scan(self, phi=np.linspace(0, 359, 359), theta=0, level=None, kinetic_energy=None, data=None): """Computes an azimuthal scan of the emitted photoelectrons. @@ -858,13 +885,13 @@ class _PED(_MSCALCULATOR): """ data = self._get_scan(scan_type='phi', level=level, theta=theta, phi=phi, kinetic_energy=kinetic_energy, data=data) - return data + return data - def get_theta_phi_scan(self, phi=np.linspace(0, 360), - theta=np.linspace(0, 90, 45), level=None, + def get_theta_phi_scan(self, phi=np.linspace(0, 360), + theta=np.linspace(0, 90, 45), level=None, kinetic_energy=None, data=None): """Computes a stereographic scan of the emitted photoelectrons. - + The azimuth ranges from 0 to 360° and the polar angle ranges from 0 to 90°. @@ -880,7 +907,7 @@ class _PED(_MSCALCULATOR): data = self._get_scan(scan_type='theta_phi', level=level, theta=theta, phi=phi, kinetic_energy=kinetic_energy, data=data, malloc={'NPH_M': 8000}) - return data + return data class _EIG(_MSCALCULATOR): @@ -946,7 +973,7 @@ class _EIG(_MSCALCULATOR): results_fname = os.path.join(self.tmp_folder, 'output/results.dat') data = self.specio.load_results(results_fname) - + energy, spec_rad = data dset.add_columns(energy=energy, spectral_radius=spec_rad) @@ -960,7 +987,7 @@ class _EIG(_MSCALCULATOR): # add the cluster to the dataset self.add_cluster_to_dset(dset) - + return self.iodata @@ -979,7 +1006,7 @@ def MSSPEC(spectroscopy='PED', **kwargs): cls = getattr(module, '_' + spectroscopy) return cls(**kwargs) - + if __name__ == "__main__": pass diff --git a/src/msspec/config.py b/src/msspec/config.py deleted file mode 100644 index eaa3047..0000000 --- a/src/msspec/config.py +++ /dev/null @@ -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 - - diff --git a/src/msspec/create_tests_results.py b/src/msspec/create_tests_results.py index d0bdd30..f0996d9 100644 --- a/src/msspec/create_tests_results.py +++ b/src/msspec/create_tests_results.py @@ -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 . +# +# Source file : src/msspec/create_tests_results.py +# Last modified: ven. 10 avril 2020 17:29:16 +# Committed by : "Sylvain Tricot " + + from msspec.tests import create_tests_results + + create_tests_results() diff --git a/src/msspec/iodata.py b/src/msspec/iodata.py index 41dbaed..d9c894d 100644 --- a/src/msspec/iodata.py +++ b/src/msspec/iodata.py @@ -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 . +# +# Source file : src/msspec/iodata.py +# Last modified: ven. 10 avril 2020 17:23:11 +# Committed by : "Sylvain Tricot " + + """ Module iodata ============= @@ -24,7 +45,7 @@ Here is an example of how to store values in a Data object: data = Data('all my data') # and append a new DataSet with its title dset = data.add_dset('Dataset 0') - + # To feed the DataSet with columns, use the add_columns method # and provide as many keywords as you like. Each key being the # column name and each value being an array holding the column @@ -47,38 +68,34 @@ Here is an example of how to store values in a Data object: import os -import numpy as np -import h5py -from lxml import etree -import msspec -from msspec.misc import LOGGER -import ase.io +import sys +from distutils.version import LooseVersion +from distutils.version import StrictVersion from io import StringIO -import wx +import ase.io +import h5py +import numpy as np import wx.grid - +from lxml import etree from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigureCanvas 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 terminaltables import AsciiTable -from distutils.version import StrictVersion, LooseVersion -import sys -#sys.path.append('../../MsSpecGui/msspecgui/msspec/gui') +import msspec 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): - # mix the values of existing theta and new theta and return the + # mix the values of existing theta and new theta and return the # unique values newx = np.linspace(np.min(x), np.max(x), nx) newy = np.linspace(np.min(y), np.max(y), ny) ux = np.unique(np.append(x, newx)) uy = np.unique(np.append(y, newy)) - + # create an empty matrix to hold the results zz = np.empty((len(ux), len(uy))) zz[:] = np.nan @@ -87,46 +104,46 @@ def cols2matrix(x, y, z, nx=88*1+1, ny=360*1+1): i = np.argwhere(ux == p[0]) j = np.argwhere(uy == p[1]) zz[i, j] = p[2] - + for i in range(len(ux)): #ok, = np.where(-np.isnan(zz[i,:])) ok, = np.where(~np.isnan(zz[i, :])) if len(ok) > 0: - xp = uy[ok] + xp = uy[ok] fp = zz[i, ok] zz[i,:] = np.interp(uy, xp, fp) - + for i in range(len(uy)): #ok, = np.where(-np.isnan(zz[:,i])) ok, = np.where(~np.isnan(zz[:, i])) if len(ok) > 0: - xp = ux[ok] + xp = ux[ok] fp = zz[ok, i] zz[:,i] = np.interp(ux, xp, fp) - + return ux, uy, zz class _DataPoint(dict): def __init__(self, *args, **kwargs): dict.__init__(self, *args, **kwargs) - + def __getattr__(self, name): if name in list(self.keys()): return self[name] else: raise AttributeError("'{}' object has no attribute '{}'".format( self.__class__.__name__, name)) - + class DataSet(object): - """ + """ This class can create an object to hold column-oriented data. :param title: The text used to entitled the dataset :type title: str :param notes: Some comments to add to the data :type notes: str - + """ def __init__(self, title, notes=""): self.title = title @@ -134,23 +151,23 @@ class DataSet(object): self._views = [] self._parameters = [] self.attributes = {} - - + + self._col_names = [] self._col_arrays = [] - self._defaults = {'bool': False, 'str': '', 'int': 0, 'float': 0., + self._defaults = {'bool': False, 'str': '', 'int': 0, 'float': 0., 'complex': complex(0)} - self._formats = {bool: '{:s}', str: '{:s}', int: '{:<20d}', + self._formats = {bool: '{:s}', str: '{:s}', int: '{:<20d}', float: '{:<20.10e}', complex: 's'} - + def _empty_array(self, val): if isinstance(val, str): t = 'S256' else: t = np.dtype(type(val)) - + if isinstance(val, bool): - default = self._defaults['bool'] + default = self._defaults['bool'] elif isinstance(val, str): default = self._defaults['str'] elif isinstance(val, int): @@ -161,9 +178,9 @@ class DataSet(object): default = self._defaults['complex'] else: raise TypeError('Not a supported type') - + return np.array([default]*len(self), dtype=t) - + def add_row(self, **kwargs): """Add a row of data into the dataset. @@ -180,17 +197,17 @@ class DataSet(object): arr = self._col_arrays[i] arr = np.append(arr, v) self._col_arrays[i] = arr - + def add_columns(self, **kwargs): """ Add columns to the dataset. - - You can provide as many columns as you want to this function. This + + You can provide as many columns as you want to this function. This function can be called several times on the same dataset but each time with different column names. Column names are given as keywords. - + :Example: - + >>> from iodata import DataSet >>> dset = DataSet('My Dataset', notes="Just an example") >>> xdata = range(10) @@ -211,7 +228,7 @@ class DataSet(object): >>> | 8 64 | >>> | 9 81 | >>> +-------+ - + """ for k, vv in list(kwargs.items()): assert k not in self._col_names, ("'{}' column already exists" @@ -228,11 +245,11 @@ class DataSet(object): def delete_rows(self, itemspec): """ Delete the rows specified with itemspec. - + """ for i in range(len(self._col_names)): self._col_arrays[i] = np.delete(self._col_arrays[i], itemspec) - + def delete_columns(self, *tags): """ Removes all columns name passed as arguments @@ -245,7 +262,7 @@ class DataSet(object): i = self._col_names.index(tag) self._col_names.pop(i) self._col_arrays.pop(i) - + def columns(self): """ Get all the column names. @@ -255,7 +272,7 @@ class DataSet(object): """ return self._col_names - + def add_view(self, name, **plotopts): """ Creates a new view named *name* with specied plot options. @@ -268,12 +285,12 @@ class DataSet(object): """ if isinstance(name, str): v = _DataSetView(self, name, **plotopts) - else: + else: v = name v.dataset = self self._views.append(v) return v - + def views(self): """Returns all the defined views in the dataset. @@ -281,7 +298,7 @@ class DataSet(object): :rtype: List of :py:class:`iodata._DataSetView` """ return self._views - + def add_parameter(self, **kwargs): """Add a parameter to store with the dataset. @@ -379,20 +396,20 @@ class DataSet(object): #fd.write(' ') fd.write(fmt.format(value)) #fd.write(str(value) + ', ') - fd.write('\n') - + fd.write('\n') + def __getitem__(self, itemspec): if isinstance(itemspec, str): return getattr(self, itemspec) title = 'untitled' new = DataSet(title) - + new._col_names = self.columns() for arr in self._col_arrays: new._col_arrays.append(np.array(arr[itemspec]).flatten()) return new - + def __setstate__(self, state): self.__dict__ = state @@ -409,7 +426,7 @@ class DataSet(object): def __iter__(self): for i in range(len(self)): - _ = {k: arr[i] for k, arr in zip(self._col_names, + _ = {k: arr[i] for k, arr in zip(self._col_names, self._col_arrays)} point = _DataPoint(_) yield point @@ -427,17 +444,17 @@ class DataSet(object): ncols = min(max_col, len(self._col_arrays)) table_data = [self._col_names[:ncols]] table_data[0].insert(0, "") - + all_indices = np.arange(0, len(self)) indices = all_indices if len(self) > max_len: indices = list(range(int(max_len/2))) + list(range(int(-max_len/2), 0)) - + _i = 0 for i in indices: if i < _i: row = ['...' for _ in range(ncols + 1)] - table_data.append(row) + table_data.append(row) row = [str(all_indices[i]),] for j in range(ncols): arr = self._col_arrays[j] @@ -446,13 +463,13 @@ class DataSet(object): row.append('...') table_data.append(row) _i = i - - table = AsciiTable(table_data) + + table = AsciiTable(table_data) table.outer_border = True table.title = self.title - table.inner_column_border = False + table.inner_column_border = False return table.table - + def __repr__(self): s = "<{}('{}')>".format(self.__class__.__name__, self.title) return s @@ -497,7 +514,7 @@ class Data(object): i = titles.index(title) self._datasets.pop(i) self._dirty = True - + def get_last_dset(self): """Get the last DataSet of the Data object. @@ -505,7 +522,7 @@ class Data(object): :rtype: :py:class:`iodata.DataSet` """ return self._datasets[-1] - + def is_dirty(self): """Wether the Data object needs to be saved. @@ -513,7 +530,7 @@ class Data(object): :rtype: bool """ return self._dirty - + def save(self, filename, append=False): """Saves the current Data to the hard drive. @@ -549,13 +566,13 @@ class Data(object): dset.title, os.path.abspath(filename))) continue grp = data_grp.create_group(dset.title) - grp.attrs['notes'] = dset.notes + grp.attrs['notes'] = dset.notes for col_name in dset.columns(): data = dset[col_name] grp.create_dataset(col_name, data=data) meta_grp.attrs['version'] = msspec.__version__ - + root = etree.Element('metainfo') # xmlize views for dset in self._datasets: @@ -563,7 +580,7 @@ class Data(object): for view in dset.views(): view_el = etree.fromstring(view.to_xml()) views_node.append(view_el) - + # xmlize parameters for dset in self._datasets: param_node = etree.SubElement(root, 'parameters', dataset=dset.title) @@ -580,7 +597,7 @@ class Data(object): meta_grp.create_dataset('info', data=np.array((xml_str,)).view('S1')) self._dirty = False LOGGER.info('Data saved in {}'.format(os.path.abspath(filename))) - + @staticmethod def load(filename): """Loads an HDF5 file from the disc. @@ -594,7 +611,7 @@ class Data(object): with h5py.File(filename, 'r') as fd: parameters = {} views = {} - + output.title = fd['DATA'].attrs['title'] for dset_name in fd['DATA'] : parameters[dset_name] = [] @@ -603,7 +620,7 @@ class Data(object): dset.notes = fd['DATA'][dset_name].attrs['notes'] for h5dset in fd['DATA'][dset_name]: dset.add_columns(**{h5dset: fd['DATA'][dset_name][h5dset].value}) - + try: vfile = LooseVersion(fd['MsSpec viewer metainfo'].attrs['version']) if vfile > LooseVersion(msspec.__version__): @@ -614,31 +631,31 @@ class Data(object): dset_name = elt0.attrib['dataset'] for elt1 in elt0.iter('parameter'): parameters[dset_name].append(elt1.attrib) - + for elt0 in root.iter('views'): dset_name = elt0.attrib['dataset'] for elt1 in elt0.iter('view'): view = _DataSetView(None, "") view.from_xml(etree.tostring(elt1)) views[dset_name].append(view) - + except Exception as err: - print(err) - - + print(err) + + for dset in output: for v in views[dset.title]: dset.add_view(v) for p in parameters[dset.title]: dset.add_parameter(**p) - output._dirty = False - return output - + output._dirty = False + return output + def __iter__(self): for dset in self._datasets: yield dset - + def __getitem__(self, key): try: titles = [d.title for d in self._datasets] @@ -653,11 +670,11 @@ class Data(object): def __str__(self): s = str([dset.title for dset in self._datasets]) return s - + def __repr__(self): s = "".format(self.title) return s - + def view(self): """Pops up a grphical window to show all the defined views of the Data object. @@ -691,8 +708,8 @@ class _DataSetView(object): legend = kwargs.get('legend', '') self._selection_conditions.append(condition) self._selection_tags.append(args) - self._plotopts['legend'].append(legend) - + self._plotopts['legend'].append(legend) + def tags(self): return self._selection_tags @@ -704,11 +721,11 @@ class _DataSetView(object): # replace all occurence of tags for tag in self.dataset.columns(): condition = condition.replace(tag, "p['{}']".format(tag)) - + for i, p in enumerate(self.dataset): if eval(condition): indices.append(i) - + values = [] for tag in tags: values.append(getattr(self.dataset[indices], tag)) @@ -724,20 +741,20 @@ class _DataSetView(object): 'plotopts': self._plotopts } root = etree.Element('root') - + return data def to_xml(self): plotopts = self._plotopts.copy() legends = plotopts.pop('legend') - + root = etree.Element('view', name=self.title) for key, value in list(plotopts.items()): root.attrib[key] = str(value) #root.attrib['dataset_name'] = self.dataset.title - - for tags, cond, legend in zip(self._selection_tags, - self._selection_conditions, + + for tags, cond, legend in zip(self._selection_tags, + self._selection_conditions, legends): curve = etree.SubElement(root, 'curve') curve.attrib['legend'] = legend @@ -745,10 +762,10 @@ class _DataSetView(object): axes = etree.SubElement(curve, 'axes') for tag in tags: variable = etree.SubElement(axes, 'axis', name=tag) - - + + return etree.tostring(root, pretty_print=False) - + def from_xml(self, xmlstr): root = etree.fromstring(xmlstr) self.title = root.attrib['name'] @@ -777,15 +794,15 @@ class _DataSetView(object): for var in curve.iter('axis'): variables.append(var.attrib['name']) tags.append(tuple(variables)) - + self._selection_conditions = conditions self._selection_tags = tags - self._plotopts['legend'] = legends - + self._plotopts['legend'] = legends + def __repr__(self): s = "<{}('{}')>".format(self.__class__.__name__, self.title) return s - + def __str__(self): try: dset_title = self.dataset.title @@ -803,7 +820,7 @@ class _GridWindow(wx.Frame): title = 'Data: ' + dset.title wx.Frame.__init__(self, parent, title=title, size=(640, 480)) self.create_grid(dset) - + def create_grid(self, dset): grid = wx.grid.Grid(self, -1) grid.CreateGrid(len(dset), len(dset.columns())) @@ -811,13 +828,13 @@ class _GridWindow(wx.Frame): grid.SetColLabelValue(ic, c) for iv, v in enumerate(dset[c]): grid.SetCellValue(iv, ic, str(v)) - + class _ParametersWindow(wx.Frame): def __init__(self, dset, parent=None): title = 'Parameters: ' + dset.title wx.Frame.__init__(self, parent, title=title, size=(400, 480)) self.create_tree(dset) - + def create_tree(self, dset): datatree = {} for p in dset.parameters(): @@ -829,10 +846,10 @@ class _ParametersWindow(wx.Frame): #group.append("{:s} = {:s}".format(p['name'], strval)) group.append("{} = {} {}".format(p['name'], p['value'], p['unit'])) datatree[p['group']] = group - + tree = wx.TreeCtrl(self, -1) root = tree.AddRoot('Parameters') - + for key in list(datatree.keys()): item0 = tree.AppendItem(root, key) for item in datatree[key]: @@ -843,7 +860,7 @@ class _ParametersWindow(wx.Frame): class _DataWindow(wx.Frame): def __init__(self, data): assert isinstance(data, (Data, DataSet)) - + if isinstance(data, DataSet): dset = data data = Data() @@ -851,9 +868,9 @@ class _DataWindow(wx.Frame): self.data = data self._filename = None self._current_dset = None - + wx.Frame.__init__(self, None, title="", size=(640, 480)) - + self.Bind(wx.EVT_CLOSE, self.on_close) # Populate the menu bar @@ -870,18 +887,18 @@ class _DataWindow(wx.Frame): sizer = wx.BoxSizer(wx.VERTICAL) #sizer.Add(self.notebook) self.SetSizer(sizer) - + self.Bind(wx.EVT_NOTEBOOK_PAGE_CHANGED, self.on_page_changed) - + self.create_notebooks() self.update_title() - + def create_notebooks(self): for key in list(self.notebooks.keys()): nb = self.notebooks.pop(key) nb.Destroy() - + for dset in self.data: nb = wx.Notebook(self, -1) self.notebooks[dset.title] = nb @@ -890,9 +907,9 @@ class _DataWindow(wx.Frame): self.create_page(nb, view) self.create_menu() - + self.show_dataset(self.data[0].title) - + def create_menu(self): menubar = wx.MenuBar() @@ -903,24 +920,24 @@ class _DataWindow(wx.Frame): menu1.Append(140, "Export\tCtrl+E") menu1.AppendSeparator() menu1.Append(199, "Close\tCtrl+Q") - + menu2 = wx.Menu() for i, dset in enumerate(self.data): menu_id = 201 + i menu2.AppendRadioItem(menu_id, dset.title) self.Bind(wx.EVT_MENU, self.on_menu_dataset, id=menu_id) - + self.Bind(wx.EVT_MENU, self.on_open, id=110) self.Bind(wx.EVT_MENU, self.on_save, id=120) self.Bind(wx.EVT_MENU, self.on_saveas, id=130) self.Bind(wx.EVT_MENU, self.on_close, id=199) - - + + menu3 = wx.Menu() menu3.Append(301, "Data") menu3.Append(302, "Cluster") menu3.Append(303, "Parameters") - + self.Bind(wx.EVT_MENU, self.on_viewdata, id=301) self.Bind(wx.EVT_MENU, self.on_viewcluster, id=302) self.Bind(wx.EVT_MENU, self.on_viewparameters, id=303) @@ -929,24 +946,24 @@ class _DataWindow(wx.Frame): menubar.Append(menu2, "&Datasets") menubar.Append(menu3, "&View") self.SetMenuBar(menubar) - + def on_open(self, event): if self.data.is_dirty(): mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do ' 'you wish to save before opening' 'another file ?'), - 'Warning: Unsaved data', + 'Warning: Unsaved data', wx.YES_NO | wx.ICON_WARNING) if mbx.ShowModal() == wx.ID_YES: self.on_saveas(wx.Event()) mbx.Destroy() - + wildcard = "HDF5 files (*.hdf5)|*.hdf5" dlg = wx.FileDialog( - self, message="Open a file...", defaultDir=os.getcwd(), + self, message="Open a file...", defaultDir=os.getcwd(), defaultFile="", wildcard=wildcard, style=wx.FD_OPEN ) - + if dlg.ShowModal() == wx.ID_OK: path = dlg.GetPath() self._filename = path @@ -954,41 +971,41 @@ class _DataWindow(wx.Frame): self.create_notebooks() dlg.Destroy() self.update_title() - + def on_save(self, event): if self._filename: if self.data.is_dirty(): self.data.save(self._filename) else: self.on_saveas(event) - + def on_saveas(self, event): overwrite = True wildcard = "HDF5 files (*.hdf5)|*.hdf5|All files (*.*)|*.*" dlg = wx.FileDialog( - self, message="Save file as ...", defaultDir=os.getcwd(), - defaultFile='{}.hdf5'.format(self.data.title.replace(' ','_')), + self, message="Save file as ...", defaultDir=os.getcwd(), + defaultFile='{}.hdf5'.format(self.data.title.replace(' ','_')), wildcard=wildcard, style=wx.FD_SAVE) dlg.SetFilterIndex(0) - + if dlg.ShowModal() == wx.ID_OK: path = dlg.GetPath() if os.path.exists(path): mbx = wx.MessageDialog(self, ('This file already exists. ' 'Do you wish to overwrite it ?'), - 'Warning: File exists', + 'Warning: File exists', wx.YES_NO | wx.ICON_WARNING) if mbx.ShowModal() == wx.ID_NO: overwrite = False mbx.Destroy() if overwrite: self.data.save(path) - self._filename = path + self._filename = path dlg.Destroy() self.update_title() - + def on_viewdata(self, event): - dset = self.data[self._current_dset] + dset = self.data[self._current_dset] frame = _GridWindow(dset, parent=self) frame.Show() @@ -1004,12 +1021,12 @@ class _DataWindow(wx.Frame): cluster_viewer.rotate_atoms(45., 45.) #cluster_viewer.show_emitter(True) win.Show() - + def on_viewparameters(self, event): dset = self.data[self._current_dset] frame = _ParametersWindow(dset, parent=self) frame.Show() - + def on_close(self, event): if self.data.is_dirty(): mbx = wx.MessageDialog(self, ('Displayed data is unsaved. Do you ' @@ -1020,15 +1037,15 @@ class _DataWindow(wx.Frame): mbx.Destroy() return self.Destroy() - + def on_menu_dataset(self, event): menu_id = event.GetId() dset_name = self.GetMenuBar().FindItemById(menu_id).GetText() self.show_dataset(dset_name) - - - def show_dataset(self, name): + + + def show_dataset(self, name): for nb in list(self.notebooks.values()): nb.Hide() self.notebooks[name].Show() @@ -1037,7 +1054,7 @@ class _DataWindow(wx.Frame): self._current_dset = name def create_page(self, nb, view): - opts = view._plotopts + opts = view._plotopts p = wx.Panel(nb, -1) figure = Figure() @@ -1060,12 +1077,12 @@ class _DataWindow(wx.Frame): toolbar.update() sizer.Add(canvas, 5, wx.ALL|wx.EXPAND) - + p.SetSizer(sizer) p.Fit() p.Show() - + for values, label in zip(view.get_data(), opts['legend']): if proj in ('ortho', 'stereo'): theta, phi, Xsec = cols2matrix(*values) @@ -1081,9 +1098,9 @@ class _DataWindow(wx.Frame): im = axes.pcolormesh(X, Y, Xsec) axes.set_yticks(R_ticks) axes.set_yticklabels(theta_ticks) - + figure.colorbar(im) - + elif proj == 'polar': values[0] = np.radians(values[0]) axes.plot(*values, label=label, picker=5, marker=opts['marker']) @@ -1107,13 +1124,13 @@ class _DataWindow(wx.Frame): axes.legend() axes.autoscale(enable=opts['autoscale']) - + # MPL events figure.canvas.mpl_connect('motion_notify_event', self.on_mpl_motion) figure.canvas.mpl_connect('pick_event', self.on_mpl_pick) - + nb.AddPage(p, view.title) - + def update_statusbar(self): sb = self.GetStatusBar() @@ -1123,7 +1140,7 @@ class _DataWindow(wx.Frame): if item.IsChecked(): sb.SetStatusText("%s" % item.GetText(), 1) break - + def update_title(self): title = "MsSpec Data Viewer" if self.data.title: @@ -1131,7 +1148,7 @@ class _DataWindow(wx.Frame): if self._filename: title += " [" + os.path.basename(self._filename) + "]" self.SetTitle(title) - + def on_mpl_motion(self, event): sb = self.GetStatusBar() try: @@ -1139,7 +1156,7 @@ class _DataWindow(wx.Frame): sb.SetStatusText(txt, 2) except Exception: pass - + def on_mpl_pick(self, event): print(event.artist) @@ -1186,7 +1203,3 @@ if __name__ == "__main__": import sys data = Data.load(sys.argv[1]) data.view() - - - - diff --git a/src/msspec/misc.py b/src/msspec/misc.py index 6054210..02786a9 100644 --- a/src/msspec/misc.py +++ b/src/msspec/misc.py @@ -1,10 +1,34 @@ -# coding: utf-8 -# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : # +#!/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 . +# +# Source file : src/msspec/misc.py +# Last modified: ven. 10 avril 2020 17:30:42 +# Committed by : "Sylvain Tricot " + + """ Module misc =========== """ + + import logging from pint import UnitRegistry import numpy as np @@ -12,6 +36,7 @@ import inspect import re import os + class XRaySource(object): MG_KALPHA = 1253.6 AL_KALPHA = 1486.6 @@ -19,10 +44,10 @@ class XRaySource(object): pass UREG = UnitRegistry() -UREG.define('rydberg = c * h * rydberg_constant = Ry') -UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0') +#UREG.define('rydberg = c * h * rydberg_constant = Ry') +#UREG.define('bohr_radius = 4 * pi * epsilon_0 * hbar**2 / electron_mass / e**2 = a0') -logging.basicConfig(level=logging.INFO) +#logging.basicConfig(level=logging.INFO) LOGGER = logging.getLogger('msspec') np.set_printoptions(formatter={'float': lambda x:'%.2f' % x}, threshold=5) diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 2bd96ed..27db224 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -1,16 +1,47 @@ -# coding: utf-8 -# vim: set et ts=4 sw=4 nu fdm=indent: +#!/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 . +# +# Source file : src/msspec/parameters.py +# Last modified: ven. 10 avril 2020 17:31:50 +# Committed by : "Sylvain Tricot " + + """ Module parameters ================= """ -import textwrap -import numpy as np + + import re -from terminaltables import AsciiTable -from msspec.misc import LOGGER, UREG, XRaySource, get_level_from_electron_configuration +import textwrap + 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): def __init__(self, name, types=None, limits=(None, None), @@ -31,7 +62,7 @@ class Parameter(object): self.fmt = fmt self.binding = binding self.private = private - self.docstring = doc + self.docstring = textwrap.dedent(doc) self._value = None self.value = default @@ -59,41 +90,41 @@ class Parameter(object): return value.to(self.unit).magnitude return value + def assert_message(self, msg, *args): + s = '\'{}\': {}'.format(self.name, msg) + s = s.format(*args) + + try: + allowed_types = '\n'.join( + ['- ' + str(_) for _ in self.allowed_types]) + except TypeError: + allowed_types = self.allowed_types + + try: + allowed_values = '\n'.join( + ['- ' + str(_) for _ in self.allowed_values]) + except TypeError: + allowed_values = self.allowed_values + + data = [ + ['PROPERTY', 'VALUE'], + ['Name', '\'{}\''.format(self.name)], + ['Allowed types', '{}'.format(allowed_types)], + ['Limits', '{} <= value <= {}'.format(self.low_limit, + self.high_limit)], + ['Unit', '{}'.format(self.unit)], + ['Allowed values', '{}'.format(allowed_values)], + ['Default', '{}'.format(self.default)] + ] + t = AsciiTable(data) + table = '\n'.join(['\t' * 2 + _ for _ in t.table.split('\n')]) + + s = "{}\n\n{}\n{}".format(s, table, self.docstring) + + return s + def check(self, value): - def assert_message(msg, *args): - s = '\'{}\': {}'.format(self.name, msg) - s = s.format(*args) - - try: - allowed_types = '\n'.join( - ['- ' + str(_) for _ in self.allowed_types]) - except TypeError: - allowed_types = self.allowed_types - - try: - allowed_values = '\n'.join( - ['- ' + str(_) for _ in self.allowed_values]) - except TypeError: - allowed_values = self.allowed_values - - data = [ - ['PROPERTY', 'VALUE'], - ['Name', '\'{}\''.format(self.name)], - ['Allowed types', '{}'.format(allowed_types)], - ['Limits', '{} <= value <= {}'.format(self.low_limit, - self.high_limit)], - ['Unit', '{}'.format(self.unit)], - ['Allowed values', '{}'.format(allowed_values)], - ['Default', '{}'.format(self.default)] - ] - t = AsciiTable(data) - table = '\n'.join(['\t' * 2 + _ for _ in t.table.split('\n')]) - - s = "{}\n\n{}\n{}".format(s, table, self.docstring) - - return s - # convert if a unit was given _value = self.convert(value) @@ -106,7 +137,7 @@ class Parameter(object): assert bool in self.allowed_types assert isinstance(_value, self.allowed_types) except AssertionError: - raise AssertionError(assert_message( + raise AssertionError(self.assert_message( 'Type {} is not an allowed type for this option ' '({} allowed)', str(type(_value)), str(self.allowed_types))) @@ -120,22 +151,22 @@ class Parameter(object): # for i, val in enumerate(_values): for val in np.array(_value).flatten(): if self.low_limit != None: - assert val >= self.low_limit, assert_message( + assert val >= self.low_limit, self.assert_message( 'Value must be >= {:s}', str(self.low_limit)) if self.high_limit != None: - assert val <= self.high_limit, assert_message( + assert val <= self.high_limit, self.assert_message( 'Value must be <= {:s}', str(self.high_limit)) - 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( + if self.allowed_values != None: # and isinstance(val, tuple(type(x) for x in self.allowed_values)): + assert val in self.allowed_values, self.assert_message( 'This value is not allowed. Please choose between ' 'one of {:s}', str(self.allowed_values)) if self.pattern != None: p = re.compile(self.pattern) m = p.match(val) - assert m != None, assert_message( + assert m != None, self.assert_message( 'Wrong string format') # _value[i] = val @@ -157,8 +188,8 @@ class Parameter(object): new_value = None try: new_value = self.binding(self) - except AttributeError: - pass + except AttributeError as err: + LOGGER.warning(err) if new_value is not None: self._value = new_value # LOGGER.debug('{:s} = {:s}'.format(self.name, str(self.value))) @@ -716,72 +747,72 @@ class GlobalParameters(BaseParameters): parameters = ( Parameter('spectroscopy', types=str, allowed_values=( 'PED', 'AED', 'LEED', 'EXAFS', 'EIG'), default='PED', - doc=textwrap.dedent(""" + doc=""" There are 4 choices for the spectroscopy option: - + - '**PED**', for Photo Electron Diffraction - '**AED**', for Auger Electron Diffraction - '**LEED**', for Low Energy Electron Diffraction - '**EXAFS**', for the Extended X-ray Absorption Fine Structure - + Additionally, a 5th keyword **EIG** is used to get deeper information about the convergence of the eigen values of multiple scattering matrix. - - The value is case insensitive. - """)), + + The value is case insensitive. + """), Parameter('algorithm', types=str, allowed_values=('expansion', 'inversion', 'correlation', 'power'), - default='expansion', doc=textwrap.dedent(""" + default='expansion', doc=""" You can choose the algorithm used for the computation of the scattering path operator. - + - '**inversion**', for the classical matrix inversion method - '**expansion**', for the Rehr-Albers series expansion - '**correlation**', for the correlation-expansion algorithm - '**power**', for the power method approximation scheme (only for spectroscopy='EIG') - + 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. - """)), + """), Parameter('polarization', types=(type(None), str), allowed_values=(None, 'linear_qOz', 'linear_xOy', - 'circular'), default=None, doc=textwrap.dedent(""" + 'circular'), default=None, doc=""" Specify the polarization of the incident light. - + - **None**, for unpolarized light - '**linear_qOz**' for a polarization vector in the :math:`(\\vec{q}0z)` plane - '**linear_xOy**' for a polarization vector in the :math:`(x0y)` plane - '**circular**' for circular dichroism - - """)), + + """), Parameter('dichroism', types=(type(None), str), 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. - """)), - Parameter('spinpol', types=bool, default=False, doc=textwrap.dedent(""" + """), + Parameter('spinpol', types=bool, default=False, doc=""" 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 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 changing it to another path. - + This folder is not automatically removed after a computation. It is removed when calling the :meth:`shutdown` calculator method: - + .. code-block:: python - + calc = MSSPEC() # the './calc' folder is created # do your calculation here calc.shutdown() # the folder is removed - - - """)) + + + """) ) BaseParameters.__init__(self) self.add_parameters(*parameters) @@ -801,7 +832,7 @@ class GlobalParameters(BaseParameters): phagen_calctype, spec_calctype = mapping[p.value] self.phagen_parameters.calctype = phagen_calctype self.spec_parameters.calctype_spectro = spec_calctype - + def bind_spinpol(self, p): if p.value == True: LOGGER.error('Spin polarization is not yet enabled in the Python version.') @@ -811,48 +842,48 @@ class GlobalParameters(BaseParameters): if p.value is not None: LOGGER.error('Dichroism is not yet enabled in the Python version.') raise NotImplemented - + class MuffintinParameters(BaseParameters): def __init__(self, phagen_parameters, spec_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 to construct the potential, is allowaed to relax around the core hole or not. - """)), - Parameter('ionicity', types=dict, default={}, doc=textwrap.dedent(""" - A dictionary to specify the ionicity of each kind of atoms. If empty, the atoms are considered to be + """), + 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 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. As an example for a LaFeO\ :sub:`3` cluster:: - + >>> calc.muffintin_parameters.ionicity = {'La': 0.15, 'Fe': 0.15, 'O': -0.1} - - 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. - - """)), + + """), Parameter('relativistic_mode', types=str, allowed_values=('non_relativistic', 'scalar_relativistic', 'spin_orbit_resolved'), - default='non_relativistic', doc=textwrap.dedent(""" + default='non_relativistic', doc=""" To tell whether to use the scalar relativistic approximation or not. - """)), + """), 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%. - """)), + """), 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 refraction at the surface. - """)), + """), Parameter('hydrogen_radius', types=(int, float), default=0.9, - limits=(0., None), unit=UREG.angstroms, doc=textwrap.dedent(""" - 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 + limits=(0., None), unit=UREG.angstroms, doc=""" + 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 radius of H atoms will be left to the program. - """)) + """) ) BaseParameters.__init__(self) self.add_parameters(*parameters) @@ -884,71 +915,71 @@ class MuffintinParameters(BaseParameters): class TMatrixParameters(BaseParameters): def __init__(self, phagen_parameters): parameters = ( - Parameter('potential', types=str, - allowed_values=('muffin_tin', 'lmto', 'spkkr', 'msf'), - default='muffin_tin', doc=textwrap.dedent(""" - This option allows to choose which kind of potential is used to compute the T-Matrix. For now, - only the internal *Muffin-Tin* potential is supported. - """)), - Parameter('potential_file', types=(str,type(None)), default=None), + Parameter('potential', types=(str, ForeignPotential), + default='muffin_tin', + doc=""" + This option allows to choose which kind of potential is + used to compute the T-Matrix. For now, only the internal + *Muffin-Tin* potential is supported.""", + ), Parameter('exchange_correlation', types=str, allowed_values=('hedin_lundqvist_real', 'hedin_lundqvist_complex', 'x_alpha_real', '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. - """)), - 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 broadening effects. - """)), + """), Parameter('lmax_mode', types=str, allowed_values=('max_ke', 'true_ke', 'imposed'), - default='true_ke', doc=textwrap.dedent(""" - This allows to control the number of basis functions used to expand the wave function on each + default='true_ke', doc=""" + This allows to control the number of basis functions used to expand the wave function on each atom. It can be: - + - '**imposed**', and will be equal to the *lmaxt* parameter (see below). It is therefore independent on both the energy and atom type. - - '**max_ke**', in this case :math:`l_{max}` is computed according to the formula + - '**max_ke**', in this case :math:`l_{max}` is computed according to the formula :math:`l_{max} = kr_{at} + 1` where :math:`k=E^{1/2}_{max}` with :math:`E_{max}` being the maximum kinetic energy. In this case :math:`l_{max}` is independent of the energy but depends on the atom number. - - '**true_ke**', in this case :math:`l_{max}` is computed according to the formula - :math:`l_{max} = kr_{at} + 1` where :math:`k=E^{1/2}_k` with :math:`E_k` being the kinetic + - '**true_ke**', in this case :math:`l_{max}` is computed according to the formula + :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. - - """)), - 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'* - """)), - 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*. For example:: - + >>> calc.tmatrix_parameters.tl_threshold = 1e-6 - + will remove all *tl* with a value :math:`< 1.10^{-6}` .. note:: - This option is compatible with any modes of the *lmax_mode* option. - - """)), - Parameter('max_tl', types=(type(None), dict), default=None, doc=textwrap.dedent(""" + This option is compatible with any modes of the *lmax_mode* option. + + """), + 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, in the case of an MgO cluster, you could write:: - + >>> calc.muffintin_parameters.max_tl = {'Mg': 20, 'O', 15} - - to tell the program to use at most 20 *tl* for Mg and 15 for O. It acts like a filter, meaning that if you - use this option, you are not required to sepcif a value for each kind of atoms in your cluster. You can + + to tell the program to use at most 20 *tl* for Mg and 15 for O. It acts like a filter, meaning that if you + use this option, you are not required to sepcif a value for each kind of atoms in your cluster. You can restrict the number of *tl* only for one type of atom for example. - + .. 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) self.add_parameters(*parameters) @@ -956,12 +987,17 @@ class TMatrixParameters(BaseParameters): self.freeze() def bind_potential(self, p): - mapping = {'muffin_tin': 'in', 'lmto': 'ex', 'spkkr': 'ex', 'msf': 'ex'} - if p.value.lower() in ('muffin_tin'): + mapping = {SPRKKRPotential: 'spkkr'} + 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' - elif p.value.lower() in ('spkkr', 'lmto', 'msf'): + elif isinstance(p.value, ForeignPotential): 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): potential = self.get_parameter('potential').value @@ -975,10 +1011,11 @@ class TMatrixParameters(BaseParameters): } self.phagen_parameters.potype = mapping[p.value] else: - self.phagen_parameters.potype = potential + mapping = {SPRKKRPotential: 'spkkr'} + self.phagen_parameters.potype = mapping[potential.value.__class__] def bind_potential_file(self, p): - pass + pass def bind_imaginery_part(self, p): self.phagen_parameters.gamma = p.value @@ -1009,37 +1046,37 @@ class SourceParameters(BaseParameters): def __init__(self, global_parameters=None, phagen_parameters=None, spec_parameters=None): parameters = ( 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. - - Common laboratories X-ray source Mg |kalpha| and Al |kalpha| lines are + + Common laboratories X-ray source Mg |kalpha| and Al |kalpha| lines are defined in the :py:class:`msspec.misc.XRaySource` class. For example: - + .. highlight:: python - + :: - - >>> from msspec.calculator import MSSPEC + + >>> from msspec.calculator import MSSPEC >>> calc = MSSPEC() >>> calc.source_parameters.energy = XRaySource.MG_KALPHA >>> print calc.source_parameters.energy 1253.6 - - - """), + + + """, default=XRaySource.MG_KALPHA), 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 - :ref:`this figure ` for questions regarding the proper + :ref:`this figure ` for questions regarding the proper orientation. - """)), + """), 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 - :ref:`this figure ` for questions regarding the proper + :ref:`this figure ` for questions regarding the proper orientation. - """)), + """), ) BaseParameters.__init__(self) self.add_parameters(*parameters) @@ -1091,13 +1128,13 @@ class DetectorParameters(BaseParameters): parameters = ( Parameter('angular_acceptance', types=(int, float), unit=UREG.degree, limits=(0., None), default=0., - doc=textwrap.dedent(""" + doc=""" The angular acceptance of the detector in degrees used when the *average_sampling* option is enabled below. - """)), + """), Parameter('average_sampling', types=(type(None), str), allowed_values=(None, 'low', 'medium', 'high'), - default=None, doc=textwrap.dedent(""" + default=None, doc=""" Used to averaged the signal over directions lying in the cone of half-angle *angular_acceptance*. The number of directions to take into account depends on the choosen @@ -1107,11 +1144,11 @@ class DetectorParameters(BaseParameters): - '**low**', to average over 5 directions - '**medium**', to average over 13 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 usual way). Otherwise, the detector is rotated. - """)) + """) ) BaseParameters.__init__(self) self.add_parameters(*parameters) @@ -1172,18 +1209,18 @@ class ScanParameters(BaseParameters): default=np.array([0.])), Parameter('kinetic_energy', types=(list, tuple, int, float), unit=UREG.eV, limits=(0., None), - default=200., doc=textwrap.dedent(""" + default=200., doc=""" if given as a list or tuple: * 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 second element as the stopping energy. - * with 3 elements, the first element is the starting energy - the second one is the stopping energy and the last one is + * with 3 elements, the first element is the starting energy + the second one is the stopping energy and the last one is the number of points. if given as a float or integer, there will be one point for the kinetic energy. - """)), + """), Parameter('ke_array', types=np.ndarray, unit=UREG.eV, default=np.array([200., ]), private=True) ) @@ -1380,58 +1417,58 @@ class CalculationParameters(BaseParameters): def __init__(self, global_parameters, phagen_parameters, spec_parameters): parameters = ( Parameter('RA_cutoff', types=int, limits=(0, 8), default=1, - doc=textwrap.dedent(""" - The Rehr-Albers cut-off parameter which controls the degree of - sphericity introduced in the description of the basis functions - used to expand the wave function around each atomic center. - It is only meaningful for the series expansion algorithm. - Its value is limited to 8 but it is rarely necessary to go beyond - 2 or 3.""")), + doc=""" + The Rehr-Albers cut-off parameter which controls the degree of + sphericity introduced in the description of the basis functions + used to expand the wave function around each atomic center. + It is only meaningful for the series expansion algorithm. + Its value is limited to 8 but it is rarely necessary to go beyond + 2 or 3."""), Parameter('scattering_order', types=int, limits=(1, 10), default=3, - doc=textwrap.dedent(""" - The scattering order. Only meaningful for the 'expansion' algorithm. - Its value is limited to 10.""")), + doc=""" + The scattering order. Only meaningful for the 'expansion' algorithm. + Its value is limited to 10."""), Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n', - 'Z_n', 'Pi_1', 'Lowdin'), + 'Z_n', 'Pi_1', 'Lowdin'), types=(type(None), str), default=None, - doc=textwrap.dedent(""" + doc=""" Enable the calculation of the coefficients for the renormalization of the multiple scattering series. - You can choose to renormalize in terms of the :math:`G_n`, the - :math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin - :math:`K^2` matrices""")), - Parameter('renormalization_omega', types=(int,float,complex), + You can choose to renormalize in terms of the :math:`G_n`, the + :math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin + :math:`K^2` matrices"""), + Parameter('renormalization_omega', types=(int,float,complex), default=1.+0j, - doc=textwrap.dedent(""" + doc=""" The :math:`\\omega` coefficient used to initialize the - renormalization alogorithm.""")), + renormalization alogorithm."""), Parameter('RA_cutoff_damping', types=int, limits=(0, None), - default=0, doc=textwrap.dedent(""" - The Rehr-Albers truncation order. If > 0, the *RA_cutoff* is - decreased by 1 every *i*\ :sup:`th` scatterer until 0, where - *i* = *RA_cutoff_damping*.""")), + default=0, doc=""" + The Rehr-Albers truncation order. If > 0, the *RA_cutoff* is + decreased by 1 every *i*\ :sup:`th` scatterer until 0, where + *i* = *RA_cutoff_damping*."""), Parameter('spin_flip', types=bool, default=False, - doc=textwrap.dedent(""" + doc=""" This parameter tells if spin-flip is authorized or not in the scattering process. - :Note: + :Note: - This option works only if the spin polarization is - enabled in your calculator object (see spinpol_).""")), + This option works only if the spin polarization is + enabled in your calculator object (see spinpol_)."""), Parameter('integrals', types=str, allowed_values=('all', 'diagonal'), - default='all', doc=textwrap.dedent(""" + default='all', doc=""" This option allows to take into account all four radial integrals (R++, R+-, R-+ and R--) in the calculation or only the diagonal radial integrals (R++ and R--) which are generally much larger. .. note:: - This option works only if the spin polarization is + This option works only if the spin polarization is enabled in your calculator object. - - """)), + + """), Parameter('path_filtering', types=(type(None), str, tuple, list), allowed_values=(None, 'forward_scattering', @@ -1439,78 +1476,78 @@ class CalculationParameters(BaseParameters): 'distance_cutoff', 'plane_wave_normal', 'plane_wave_spin_averaged'), - default=None, doc=textwrap.dedent(""" - Used to activate some filters. It is possible to specify several + default=None, doc=""" + Used to activate some filters. It is possible to specify several of them by grouping them in a tuple or a list. For example:: >>> my_filters = ('forward_scattering', 'backward_scattering') >>> calc.calculation_parameters.path_filtering = my_filters - - """)), + + """), Parameter('off_cone_events', types=int, limits=(0, None), default=1, - doc=textwrap.dedent(""" - Used in conjunction with the '*forward_scattering*' filter. - If the number of scattering processes outside the forward (or - backward) scattering cone is greater than this number, then the - path is rejected and its contribution to the scattering path - operator won’t be computed. - """)), + doc=""" + Used in conjunction with the '*forward_scattering*' filter. + If the number of scattering processes outside the forward (or + backward) scattering cone is greater than this number, then the + path is rejected and its contribution to the scattering path + operator won’t be computed. + """), Parameter('scattering_order_cutoff', types=int, limits=(0, 10), - default=2, doc=textwrap.dedent(""" - Used in conjunction with the ‘*plane_wave_normal*’ filter. It states - to activate the plane wave approximation (which is fast but - less accurate) to compute the contribution when the scattering order - is greater than this value.""")), + default=2, doc=""" + Used in conjunction with the ‘*plane_wave_normal*’ filter. It states + to activate the plane wave approximation (which is fast but + less accurate) to compute the contribution when the scattering order + is greater than this value."""), Parameter('distance', types=(int, float), limits=(0, None), - unit=UREG.angstroms, default=10., doc=textwrap.dedent(""" - Used with the '*distance_cut_off*' filter. Paths whose length is - larger than this value are simply rejected.""")), + unit=UREG.angstroms, default=10., doc=""" + Used with the '*distance_cut_off*' filter. Paths whose length is + larger than this value are simply rejected."""), Parameter('vibrational_damping', types=(type(None), str), 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: - '**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), - 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 - """)), + """), Parameter('debye_temperature', types=(int, float), limits=(0, None), - unit=UREG.degK, default=420., doc=textwrap.dedent(""" - The Debye temperature used for the calculation of the mean square - displacements if *use_debye_model* = True""")), + unit=UREG.degK, default=420., doc=""" + The Debye temperature used for the calculation of the mean square + displacements if *use_debye_model* = True"""), Parameter('use_debye_model', types=bool, default=False, - doc=textwrap.dedent(""" - No matter the way you compute the effect of atomic vibrations, - you need the mean square displacements of atoms. It can be computed + doc=""" + No matter the way you compute the effect of atomic vibrations, + you need the mean square displacements of atoms. It can be computed internally following the Debye model if you set this option to True. - """)), + """), Parameter('vibration_scaling', types=(int, float), - limits=(0., None), default=1.2, doc=textwrap.dedent(""" - Used to simulate the fact that surface atoms vibrate more than - bulk ones. It is a factor applied to the mean square displacements. - """)), + limits=(0., None), default=1.2, doc=""" + Used to simulate the fact that surface atoms vibrate more than + bulk ones. It is a factor applied to the mean square displacements. + """), Parameter('basis_functions', types=str, allowed_values=( 'plane_wave', 'spherical'), default='spherical', private=True), Parameter('cutoff_factor', types=(int, float), limits=(1e-4, 999.9999), default=0.01, private=True), Parameter('mean_free_path', types=(int, float, str), default='SeahDench', allowed_values=('mono', 'SeahDench'), - doc=textwrap.dedent(""" + doc=""" 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 the keyword 'mono' to use a formula valid only for monoelemental samples - Enter the keyword 'SeahDench' to use the Seah and Dench formula. - + .. note:: - + The mean free path is only taken into account when the input T-matrix corresponds - 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. - - """)), + + """), ) BaseParameters.__init__(self) @@ -1537,11 +1574,11 @@ class CalculationParameters(BaseParameters): 'Pi_1' : 4, 'Lowdin' : 5} # Check that the method is neither 'Z_n' nor 'K^2' for other - # 'spetroscopy' than EIG + # 'spetroscopy' than EIG try: if (self.global_parameters.spectroscopy == 'PED' and self.global_parameters.algorithm == 'expansion'): - #assert( p.value in (None, 'Sigma_n', 'G_n') ) + #assert( p.value in (None, 'Sigma_n', 'G_n') ) assert( p.value in p.allowed_values) elif (self.global_parameters.spectroscopy == 'EIG' and self.global_parameters.algorithm == 'inversion'): @@ -1561,7 +1598,7 @@ class CalculationParameters(BaseParameters): self.spec_parameters.calc_renr = omega.real self.spec_parameters.calc_reni = omega.imag LOGGER.info("Renormalization omega set to \'{}\'".format(p.value)) - + def bind_RA_cutoff_damping(self, p): self.spec_parameters.calc_ino = p.value LOGGER.info('Rehr-Albers cutoff damping set to %s', p.value) @@ -1728,8 +1765,8 @@ class PEDParameters(BaseParameters): def __init__(self, phagen_parameters, spec_parameters): parameters = ( Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$', - default='1s', doc=textwrap.dedent(""" - The level is the electronic level where the electron comes from. + default='1s', doc=""" + The level is the electronic level where the electron comes from. It is written: *nlJ* where: @@ -1742,7 +1779,7 @@ class PEDParameters(BaseParameters): >>> calc.spectroscopy_parameters.level = '2p3/2' >>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2' - """)), + """), Parameter('final_state', types=int, limits=(-1, 2), default=2), Parameter('spin_orbit', types=(type(None), str), allowed_values=(None, 'single', 'both'), default=None), @@ -1779,8 +1816,8 @@ class EIGParameters(BaseParameters): def __init__(self, phagen_parameters, spec_parameters): parameters = ( Parameter('level', types=str, pattern=r'\d+[spdfgSPDFG](\d/2)?$', - default='1s', doc=textwrap.dedent(""" - The level is the electronic level where the electron comes from. + default='1s', doc=""" + The level is the electronic level where the electron comes from. It is written: *nlJ* where: @@ -1793,21 +1830,21 @@ class EIGParameters(BaseParameters): >>> calc.spectroscopy_parameters.level = '2p3/2' >>> calc.spectroscopy_parameters.level = '2p' # is equivalent to '2p1/2' - """)), + """), Parameter('final_state', types=int, limits=(-1, 2), default=2), Parameter('spin_orbit', types=(type(None), str), 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. - """)), + """), Parameter('method', types=(str,), default='EPSI', allowed_values=['AITK', 'RICH', 'SALZ', 'EPSI', 'EPSG', 'RHOA', 'THET', 'LEGE', 'CHEB', 'OVER', 'DURB', 'DLEV', 'TLEV', 'ULEV', 'VLEV', 'ELEV', 'EULE', 'GBWT', 'VARI', 'ITHE', 'EALG'], - doc=textwrap.dedent("""The convergence acceleration scheme to be used. - """)), + doc="""The convergence acceleration scheme to be used. + """), ) BaseParameters.__init__(self) self.add_parameters(*parameters) @@ -1832,5 +1869,3 @@ class EIGParameters(BaseParameters): def bind_kernel_matrix_spectrum(self, p): value = int(p.value) self.spec_parameters.eigval_ispectrum_ne = value - - diff --git a/src/msspec/tests.py b/src/msspec/tests.py index 9cdf7c6..50246ae 100644 --- a/src/msspec/tests.py +++ b/src/msspec/tests.py @@ -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 . +# +# Source file : src/msspec/tests.py +# Last modified: ven. 10 avril 2020 17:33:28 +# Committed by : "Sylvain Tricot " -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 pickle import dumps -import unittest -import sys -import os +from ase.build import bulk + +from msspec.calculator import MSSPEC +from msspec.misc import set_log_level +from msspec.utils import * set_log_level('error') RESULTS_FILENAME = os.path.join(os.path.dirname(__file__), 'results.txt') + def perform_test(obj, funcname, filename=RESULTS_FILENAME): f = getattr(obj, '_' + funcname) output = md5(dumps(f())).hexdigest() diff --git a/src/msspec/utils.py b/src/msspec/utils.py index 18574bf..3097dd3 100644 --- a/src/msspec/utils.py +++ b/src/msspec/utils.py @@ -1,5 +1,26 @@ -# -*- encoding: utf-8 -*- -# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : # +#!/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 . +# +# Source file : src/msspec/utils.py +# Last modified: ven. 10 avril 2020 15:49:35 +# Committed by : "Sylvain Tricot " + """ Module utils @@ -8,20 +29,159 @@ Module utils """ + +import re + import numpy as np -from ase import Atoms, Atom -from ase.visualize import view +from ase import Atom +from ase import Atoms -class MsSpecAtoms(Atoms): - def __init__(self, *args, **kwargs): - Atoms.__init__(self, *args, **kwargs) - self.__absorber_index = None +from msspec.iodata import Data +from msspec.misc import LOGGER - def set_absorber(self, index): - self.__absorber_index = index - def get_absorber(self): - return self.__absorber_index +class ForeignPotential(object): + def __init__(self): + self.data = Data(title='Foreign Potential') + # Load exported potentials + # phagen_data is a dictionary with: + # self.phagen_data = { + # 'VintTotal' : , + # 'VintCoulomb': , + # 'RHOint' : , + # 'types': [ + # { + # 'Z' : , + # 'RWS' : , + # 'data': + # }, + # ... + # ... + # ... + # ] + # } + 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*(?PIQ.*)\n' + r'(?P(.*\n)+?))\*+'), + [int] + [float] * 3) + self.types_data = read(content, + (r'^\s*TYPES\s*\n(\s*(?PIT.*)\n' + r'(?P(.*\n)+?))\*+'), + [int, str] + [int] * 4) + self.occ_data = read(content, + (r'^\s*OCCUPATION\s*\n(\s*(?PIQ.*)\n' + r'(?P(.*\n)+?))\*+'), + [int] * 5 + [float]) + + for f in exported_files: + # get the IT from the filename + # m=re.match('SPRKKR-IT_(?P\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.*?)\s+' + r'VMTZ_Coulomb\s*=\s*(?P.*?)\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): def __init__(self, *args, **kwargs): @@ -30,137 +190,119 @@ class EmptySphere(Atom): def get_atom_index(atoms, x, y, z): - """ Return the index of the atom that is the closest to the coordiantes - given as parameters. + """ Return the index of the atom that is the closest to the coordiantes + given as parameters. - :param ase.Atoms atoms: an ASE Atoms object - :param float x: the x position in angstroms - :param float y: the y position in angstroms - :param float z: the z position in angstroms + :param ase.Atoms atoms: an ASE Atoms object + :param float x: the x position in angstroms + :param float y: the y position in angstroms + :param float z: the z position in angstroms - :return: the index of the atom as an integer - :rtype: int - """ - # get all distances - d = np.linalg.norm(atoms.get_positions() - np.array([x, y, z]), axis = 1) - # get the index of the min distance - i = np.argmin(d) - # return the index and the corresponding distance - return i + :return: the index of the atom as an integer + :rtype: int + """ + # get all distances + d = np.linalg.norm(atoms.get_positions() - np.array([x, y, z]), axis=1) + # get the index of the min distance + i = np.argmin(d) + # return the index and the corresponding distance + return i def center_cluster(atoms, invert=False): - """ Centers an Atoms object by translating it so the origin is roughly - at the center of the cluster. - The function supposes that the cluster is wrapped inside the unit 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 - origin is at the center of the cluster - - :param ase.Atoms atoms: an ASE Atoms object - :param bool invert: if True, performs the opposite translation (uncentering the cluster) - - """ - for i, cell_vector in enumerate(atoms.get_cell()): - if invert: - atoms.translate(0.5*cell_vector) - else: - atoms.translate(-0.5*cell_vector) + """ Centers an Atoms object by translating it so the origin is roughly + at the center of the cluster. + The function supposes that the cluster is wrapped inside the unit 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 + origin is at the center of the cluster + :param ase.Atoms atoms: an ASE Atoms object + :param bool invert: if True, performs the opposite translation + (uncentering the cluster) + """ + for i, cell_vector in enumerate(atoms.get_cell()): + if invert: + atoms.translate(0.5*cell_vector) + else: + atoms.translate(-0.5*cell_vector) def cut_sphere(atoms, radius, center=(0, 0, 0)): + """ Removes all the atoms of an Atoms object outside a sphere with a + given radius + + :param ase.Atoms atoms: an ASE Atoms object + :param float radius: the radius of the sphere + + :return: The modified atom cluster + :rtype: ase.Atoms + """ 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 - given radius - - :param ase.Atoms atoms: an ASE Atoms object - :param float radius: the radius of the sphere - - :return: The modified atom cluster - :rtype: ase.Atoms - """ - if radius is None: - raise ValueError("radius not set") - - 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): - """ Removes all the atoms of an Atoms object outside a cylinder with a - given axis and radius + """ Removes all the atoms of an Atoms object outside a cylinder with a + given axis and radius - :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 float radius: the radius of the cylinder + :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 float radius: the radius of the cylinder - :return: The modified atom cluster - :rtype: ase.Atoms - """ - if radius is None: - raise ValueError("radius not set") + :return: The modified atom cluster + :rtype: ase.Atoms + """ + if radius is None: + raise ValueError("radius not set") - new_atoms = atoms.copy() + new_atoms = atoms.copy() - dims = {"x": 0, "y": 1, "z": 2} - if axis in dims: - axis = dims[axis] - else: - raise ValueError("axis not valid, must be 'x','y', or 'z'") + dims = {"x": 0, "y": 1, "z": 2} + if axis in dims: + axis = dims[axis] + else: + raise ValueError("axis not valid, must be 'x','y', or 'z'") - del_list = [] - for index, position in enumerate(new_atoms.positions): - # calculating the distance of the atom to the given axis - r = 0 - for dim in range(3): - if dim != axis: - r = r + position[dim]**2 - r = np.sqrt(r) + del_list = [] + for index, position in enumerate(new_atoms.positions): + # calculating the distance of the atom to the given axis + r = 0 + for dim in range(3): + if dim != axis: + r = r + position[dim]**2 + r = np.sqrt(r) - if r > radius: - del_list.append(index) + if r > radius: + del_list.append(index) - del_list.reverse() - for index in del_list: - del new_atoms[index] + del_list.reverse() + for index in del_list: + del new_atoms[index] - return new_atoms - -def cut_cone(atoms, radius, z = 0): + return new_atoms + + +def cut_cone(atoms, radius, z=0): """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. :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 - :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 :return: A new cluster. :rtype: :py:class:`ase.Atoms` """ new_atoms = atoms.copy() - origin = np.array((0, 0, 0)) max_theta = np.arctan(radius/(-z)) - + u = np.array((0, 0, -z)) normu = np.linalg.norm(u) new_atoms.translate(u) @@ -168,20 +310,20 @@ def cut_cone(atoms, radius, z = 0): for i in range(len(new_atoms)): v = new_atoms[i].position normv = np.linalg.norm(v) - + _ = np.dot(u, v)/normu/normv if _ == 0: print(v) theta = np.arccos(_) if theta <= max_theta: indices.append(i) - - + new_atoms = new_atoms[indices] - new_atoms.translate(-u) # pylint: disable=invalid-unary-operand-type - + new_atoms.translate(-u) + return new_atoms + def cut_plane(atoms, x=None, y=None, z=None): """ Removes the atoms whose coordinates are higher (or lower, for a 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 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 - #y-coordinate is higher than 3.6, and/or z-coordinate is higher than 0 - #is deleted. + # every atom whose x-coordinate is higher than 5 or lower than -5, + # and/or y-coordinate is higher than 3.6, and/or z-coordinate is higher + # than 0 is deleted. :param ase.Atoms atoms: an ASE Atoms object :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_values = [x, y, z] 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)): - 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" else: try: @@ -216,7 +360,7 @@ def cut_plane(atoms, x=None, y=None, z=None): dim_values[i] = [-np.inf, value] else: dim_values[i] = [value, np.inf] - except: + except Exception: dim_values[i] = [value, value] 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) 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] 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'): - """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 :type cluster: Atoms object @@ -245,7 +392,8 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, :type diameter: float :param planes: the number of planes of your cluster :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 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 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) if (n < nmin): 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 c = cell[:, 2].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: - 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: - l = diameter + min_diameter = diameter - rep = int(2*l/min(a,c)) # number of repetition in each direction - cluster = cluster.repeat((rep, rep, rep)) # repeat the cluster + # number of repetition in each direction + rep = int(3*min_diameter/min(a, c)) - center_cluster(cluster) # center the cluster - cluster.set_cell(cell) # reset the cell - cluster = cut_plane(cluster, z=eps) # cut the cluster so that we have a centered surface + # repeat the cluster + cluster = cluster.repeat((rep, rep, rep)) - i = np.where(cluster.get_tags() == emitter_tag) # positions where atoms have the tag of the emitter_tag - 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 - all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # an array of all unique z + # center the cluster + center_cluster(cluster) - 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 - ze = all_ze.max() # the height of the emitter's plane + # reset the cell + 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 + i = np.where(cluster.get_tags() == emitter_tag) + + # 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() + + # 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] + n = np.where(all_z == all_z.max())[0][0] - np.where( + all_z == all_ze.max())[0][0] ze = all_ze.max() + # values of x and y of the emitter + tx, ty = get_xypos(cluster, ze, symbol) - tx, ty = get_xypos(cluster, ze, symbol) # values of x and y of the emitter - Atoms.translate(cluster, [-tx, -ty, 0]) # center the cluster on the emitter + # center the cluster on the emitter + Atoms.translate(cluster, [-tx, -ty, 0]) - 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 - Atoms.translate(cluster, [0, 0, -z_cut]) # translate the surface at z=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 - if planes!=0: - all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) # an array of all unique remaining z + if planes != 0: + # an array of all unique remaining z + all_z = np.sort(np.unique(np.round(cluster.get_positions()[:, 2], 4))) + zplan = all_z[-planes] xplan, yplan = get_xypos(cluster, zplan) 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) 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'): - 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: raise NameError('Unkknown shape specifier: \"{}\"'.format(shape)) - 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 - cluster = cut_plane(cluster, z=zcut) # cut the right number of planes + if planes != 0: + # calculate where to cut to get 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 - assert emitter_plane < np.alen(all_z), "There are not enough existing plans." + # cut the right number of planes + 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 Atoms.translate(cluster, [0, 0, -ze]) # put the emitter in (0,0,0) diff --git a/src/msspec/version.py b/src/msspec/version.py index 0ef25a6..ca1e1fb 100644 --- a/src/msspec/version.py +++ b/src/msspec/version.py @@ -1,10 +1,33 @@ -# coding: utf-8 -# vim: set et sw=4 ts=4 sts=4 nu ai cc=+0 fdm=indent mouse=a: +#!/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 . +# +# Source file : src/msspec/version.py +# Last modified: ven. 10 avril 2020 17:34:38 +# Committed by : "Sylvain Tricot " + -from setuptools_scm import get_version -from pkg_resources import parse_version, DistributionNotFound, get_distribution 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 # 1- Try to read it from the git info # 2- If it fails, try to read it from the distribution file diff --git a/src/setup.py b/src/setup.py index 075d91c..ac2b26a 100644 --- a/src/setup.py +++ b/src/setup.py @@ -1,22 +1,31 @@ -# coding: utf-8 -# vim: set fdm=indent ts=4 sw=4 sts=4 et tw=80 ai cc=+0 mouse=a nu : +#!/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. -from setuptools import setup, Extension, find_packages -from distutils.file_util import copy_file +# 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. -from distutils.command.build import build as _build -from setuptools.command.build_ext import build_ext as _build_ext -from setuptools.command.install import install as _install -from distutils.command.clean import clean as _clean +# You should have received a copy of the GNU General Public License +# along with this msspec. If not, see . +# +# Source file : src/setup.py +# Last modified: mar. 07 avril 2020 17:01:42 +# Committed by : "Sylvain Tricot " -import subprocess -import traceback -import os +from setuptools import setup, find_packages +from version import __version__ import sys -import glob sys.path.insert(0, "msspec") -from version import __version__ with open('setup_requirements.txt', 'r') as fd: SETUP_REQUIREMENTS = fd.read().strip().split('\n') @@ -26,46 +35,46 @@ with open('requirements.txt', 'r') as fd: if __name__ == "__main__": setup(name='msspec', - version=__version__, - include_package_data=True, - packages=find_packages(include='msspec.*'), - setup_requires=SETUP_REQUIREMENTS, - install_requires=REQUIREMENTS, + version=__version__, + include_package_data=True, + packages=find_packages(include='msspec.*'), + setup_requires=SETUP_REQUIREMENTS, + install_requires=REQUIREMENTS, - author='Didier Sébilleau, Sylvain Tricot', - author_email='sylvain.tricot@univ-rennes1.fr', - maintainer='Sylvain Tricot', - maintainer_email='sylvain.tricot@univ-rennes1.fr', - url='https://msspec.cnrs.fr', - description=('A multiple scattering package for sepectroscopies using ' - 'electrons to probe materials'), - long_description="""MsSpec is a Fortran package to compute the - cross-section of several spectroscopies involving one (or more) - electron(s) as the probe. This package provides a python interface to - control all the steps of the calculation. - - Available spectroscopies: - * Photoelectron diffraction - * Auger electron diffraction - * Low energy electron diffraction - * X-Ray absorption spectroscopy - * Auger Photoelectron coincidence spectroscopy - * Computation of the spectral radius""", - download_url='https://msspec.cnrs.fr/downloads.html', - # See https://pypi.python.org/pypi?%3Aaction=list_classifiers - classifiers=[ - 'Development Status :: 3 - Alpha', - 'Environment :: Console', - 'Intended Audience :: Science/Research', - 'License :: OSI Approved :: GNU General Public License (GPL)', - 'Natural Language :: English', - 'Operating System :: Microsoft :: Windows :: Windows 10', - 'Operating System :: POSIX :: Linux', - 'Operating System :: MacOS :: MacOS X', - 'Programming Language :: Fortran', - 'Programming Language :: Python :: 3 :: Only', - 'Topic :: Scientific/Engineering :: Physics', + author='Didier Sébilleau, Sylvain Tricot', + author_email='sylvain.tricot@univ-rennes1.fr', + maintainer='Sylvain Tricot', + maintainer_email='sylvain.tricot@univ-rennes1.fr', + url='https://msspec.cnrs.fr', + description=('A multiple scattering package for sepectroscopies ' + 'using electrons to probe materials'), + long_description="""MsSpec is a Fortran package to compute the + cross-section of several spectroscopies involving one (or more) + electron(s) as the probe. This package provides a python interface to + control all the steps of the calculation. + + Available spectroscopies: + * Photoelectron diffraction + * Auger electron diffraction + * Low energy electron diffraction + * X-Ray absorption spectroscopy + * Auger Photoelectron coincidence spectroscopy + * Computation of the spectral radius""", + download_url='https://msspec.cnrs.fr/downloads.html', + # See https://pypi.python.org/pypi?%3Aaction=list_classifiers + classifiers=[ + 'Development Status :: 3 - Alpha', + 'Environment :: Console', + 'Intended Audience :: Science/Research', + 'License :: OSI Approved :: GNU General Public License (GPL)', + 'Natural Language :: English', + 'Operating System :: Microsoft :: Windows :: Windows 10', + 'Operating System :: POSIX :: Linux', + 'Operating System :: MacOS :: MacOS X', + 'Programming Language :: Fortran', + 'Programming Language :: Python :: 3 :: Only', + 'Topic :: Scientific/Engineering :: Physics', ], - keywords='spectroscopy atom electron photon multiple scattering', - license='GPL', + keywords='spectroscopy atom electron photon multiple scattering', + license='GPL', )