amd._reader._Reader._construct_periodic_set() - Code Metrics - Inspection of "io clean 3" - dwiddo/average-minimum-distance - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( ca5a49...2d0e3f )

by Daniel

created 2022-04-15 17:00 UTC

amd._reader._Reader._construct_periodic_set() A

↳ Parent: amd._reader

Complexity

Conditions

Size

Total Lines	18
Code Lines	13

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	13
dl	0
loc	18
rs	9.75
c	0
b	0
f	0
cc	2
nop	6

"""Contains base reader class for the io module. 

This class implements the converters from CifBlock, Entry to PeriodicSets.
"""

import warnings
from typing import Callable, Iterable, Sequence, Tuple

import numpy as np
import ase.spacegroup.spacegroup    # parse_sitesym
import ase.io.cif

from .periodicset import PeriodicSet
from .utils import cellpar_to_cell


class _Reader:

    """Base Reader class. Contains parsers for converting ase CifBlock
    and ccdc Entry objects to PeriodicSets.

    Intended use:

    First make a new method for _Reader converting object to PeriodicSet
    (e.g. named _X_to_PSet). Then make this class outline:

    class XReader(_Reader):
        def __init__(self, ..., **kwargs):

        super().__init__(**kwargs)

        # setup and checks

        # make 'iterable' which yields objects to be converted (e.g. CIFBlock, Entry)

        # set self._generator like this
        self._generator = self._read(iterable, self._X_to_PSet)
    """

    disorder_options = {'skip', 'ordered_sites', 'all_sites'}
    reserved_tags = {
        'motif',
        'cell',
        'name',
        'asymmetric_unit',
        'wyckoff_multiplicities',
        'types',
        'filename',}
    atom_site_fract_tags = [
        '_atom_site_fract_x',
        '_atom_site_fract_y',
        '_atom_site_fract_z',]
    atom_site_cartn_tags = [
        '_atom_site_cartn_x',
        '_atom_site_cartn_y',
        '_atom_site_cartn_z',]
    symop_tags = [
        '_space_group_symop_operation_xyz',
        '_space_group_symop.operation_xyz',
        '_symmetry_equiv_pos_as_xyz',]

    equiv_site_tol = 1e-3

    def __init__(

            self,
            remove_hydrogens=False,
            disorder='skip',
            heaviest_component=False,
            show_warnings=True,
            extract_data=None,
            include_if=None):

        if disorder not in _Reader.disorder_options:
            raise ValueError(f'disorder parameter {disorder} must be one of {_Reader.disorder_options}')


        if extract_data is None:
            self.extract_data = {}
        else:
            _validate_extract_data(extract_data)
            self.extract_data = extract_data

        if include_if is None:
            self.include_if = ()
        elif not all(callable(func) for func in include_if):
            raise ValueError('include_if must be a list of callables')
        else:
            self.include_if = include_if
        

        self.remove_hydrogens = remove_hydrogens
        self.disorder = disorder
        self.heaviest_component = heaviest_component
        self.show_warnings = show_warnings
        self.current_name = None
        self.current_filename = None
        self._generator = []

    def __iter__(self):
        yield from self._generator

    def read_one(self):
        """Read the next (or first) item."""
        return next(iter(self._generator))

    def _map(self, func: Callable, iterable: Iterable) -> Iterable[PeriodicSet]:
        """Iterates over iterable, passing items through parser and yielding the 

        result if it is not None. Applies warning and include_if filter.
        """

        with warnings.catch_warnings():
            if not self.show_warnings:
                warnings.simplefilter('ignore')
            for item in iterable:
                if all(check(item) for check in self.include_if):
                    res = func(item)
                    if res is not None:
                        yield res

    def _cifblock_to_periodicset(self, block) -> PeriodicSet:
        """ase.io.cif.CIFBlock --> PeriodicSet. Returns None for a "bad" set."""

        self.current_name = block.name
        asym_unit, asym_symbols, sitesym, cell = self._cifblock_to_asym_unit(block)

        # indices of sites to remove
        remove = []
        if self.remove_hydrogens:
            remove.extend((i for i, sym in enumerate(asym_symbols) if sym in 'HD'))


        # find disordered sites
        asym_is_disordered = []
        occupancies = block.get('_atom_site_occupancy')
        labels = block.get('_atom_site_label')
        if occupancies is not None:
            disordered = []     # indices where there is disorder
            for i, (occ, label) in enumerate(zip(occupancies, labels)):
                if _atom_has_disorder(label, occ):
                    if i not in remove:
                        disordered.append(i)
                        asym_is_disordered.append(True)
                else:
                    asym_is_disordered.append(False)

            if self.disorder == 'skip' and len(disordered) > 0:
                warnings.warn(f'Skipping {self.current_name} as structure is disordered')
                return None

            if self.disorder == 'ordered_sites':
                remove.extend(disordered)

        # remove sites
        asym_unit = np.delete(asym_unit, remove, axis=0)
        asym_symbols = [s for i, s in enumerate(asym_symbols) if i not in remove]
        asym_is_disordered = [v for i, v in enumerate(asym_is_disordered) if i not in remove]

        keep_sites = self._validate_sites(asym_unit, asym_is_disordered)
        if keep_sites is not None:
            asym_unit = asym_unit[keep_sites]
            asym_symbols = [sym for sym, keep in zip(asym_symbols, keep_sites) if keep]

        if self._has_no_valid_sites(asym_unit):
            return None

        data = {key: func(block) for key, func in self.extract_data.items()}
        periodic_set = self._construct_periodic_set(asym_unit, asym_symbols, sitesym, cell, **data)
        return periodic_set

    def _cifblock_to_asym_unit(self, block):
        """ase.io.cif.CIFBlock --> 

        asymmetric unit (frac coords), asym_symbols, cell, symops (as strings)"""
        

        cell = block.get_cell().array
        asym_unit = [block.get(name) for name in _Reader.atom_site_fract_tags]
        if None in asym_unit:
            asym_motif = [block.get(name) for name in _Reader.atom_site_cartn_tags]
            if None in asym_motif:
                warnings.warn(f'Skipping {self.current_name} as coordinates were not found')
                return None
            asym_unit = np.array(asym_motif) @ np.linalg.inv(cell)
        asym_unit = np.mod(np.array(asym_unit).T, 1)
        

        try:
            asym_symbols = block.get_symbols()
        except ase.io.cif.NoStructureData as _:
            asym_symbols = ['Unknown' for _ in range(len(asym_unit))]
        

        sitesym = ['x,y,z', ]
        for tag in _Reader.symop_tags:
            if tag in block:
                sitesym = block[tag]
                break

        if isinstance(sitesym, str):
            sitesym = [sitesym]
        

        return asym_unit, asym_symbols, sitesym, cell

    def _entry_to_periodicset(self, entry) -> PeriodicSet:
        """ccdc.entry.Entry --> PeriodicSet. Returns None for a "bad" set."""

        self.current_name = entry.identifier
        crystal = entry.crystal

        if not entry.has_3d_structure:
            warnings.warn(f'Skipping {self.current_name} as entry has no 3D structure')
            return None

        # first disorder check, if skipping. If occ == 1 for all atoms but the entry
        # or crystal is listed as having disorder, skip (can't know where disorder is).
        # If occ != 1 for any atoms, we wait to see if we remove them before skipping.
        molecule = crystal.disordered_molecule
        if self.disorder == 'ordered_sites':
            molecule.remove_atoms(a for a in molecule.atoms if a.label.endswith('?'))

        may_have_disorder = False
        if self.disorder == 'skip':
            for a in molecule.atoms:
                occ = a.occupancy
                if _atom_has_disorder(a.label, occ):
                    may_have_disorder = True
                    break

            if not may_have_disorder:
                if crystal.has_disorder or entry.has_disorder:
                    warnings.warn(f'Skipping {self.current_name} as structure is disordered')
                    return None

        # make same as cifblock version??
        if self.remove_hydrogens:
            molecule.remove_atoms(a for a in molecule.atoms if a.atomic_symbol in 'HD')

        if self.heaviest_component and len(molecule.components) > 1:
            molecule = _heaviest_component(molecule)

        crystal.molecule = molecule

        # by here all atoms to be removed have been (except via ordered_sites).
        # If disorder == 'skip' and there were atom(s) with occ < 1 found
        # eariler, we check if all such atoms were removed. If not, skip.
        if self.disorder == 'skip' and may_have_disorder:
            for a in crystal.disordered_molecule.atoms:
                occ = a.occupancy
                if _atom_has_disorder(a.label, occ):
                    warnings.warn(f'Skipping {self.current_name} as structure is disordered')
                    return None

        # if disorder is all_sites, we need to know where disorder is to ignore overlaps
        asym_is_disordered = []     # True/False list same length as asym unit
        if self.disorder == 'all_sites':
            for a in crystal.asymmetric_unit_molecule.atoms:
                occ = a.occupancy
                if _atom_has_disorder(a.label, occ):
                    asym_is_disordered.append(True)
                else:
                    asym_is_disordered.append(False)

        # check all atoms have coords. option/default remove unknown sites?
        if not molecule.all_atoms_have_sites or \
           any(a.fractional_coordinates is None for a in molecule.atoms):
            warnings.warn(f'Skipping {self.current_name} as some atoms do not have sites')
            return None

        asym_unit, asym_symbols, sitesym, cell = self._crystal_to_asym_unit(crystal)

        # remove overlapping sites, check sites exist
        keep_sites = self._validate_sites(asym_unit, asym_is_disordered)
        if keep_sites is not None:
            asym_unit = asym_unit[keep_sites]
            asym_symbols = [sym for sym, keep in zip(asym_symbols, keep_sites) if keep]

        if self._has_no_valid_sites(asym_unit):
            return None

        entry.crystal.molecule = crystal.disordered_molecule
        data = {key: func(entry) for key, func in self.extract_data.items()}
        periodic_set = self._construct_periodic_set(asym_unit, asym_symbols, sitesym, cell, **data)
        return periodic_set

    def _crystal_to_asym_unit(self, crystal):
class Foo:
    def some_method(self, x, y):
        return x + y;
        """ase.io.cif.CIFBlock -->
        asymmetric unit (frac coords), asym_symbols, symops, cell"""
        

        asym_atoms = crystal.asymmetric_unit_molecule.atoms
        asym_unit = np.array([tuple(a.fractional_coordinates) for a in asym_atoms])
        asym_unit = np.mod(asym_unit, 1)
        asym_symbols = [a.atomic_symbol for a in asym_atoms]
        cell = cellpar_to_cell(*crystal.cell_lengths, *crystal.cell_angles)
        sitesym = crystal.symmetry_operators
        if not sitesym:
            sitesym = ['x,y,z', ]
        return asym_unit, asym_symbols, sitesym, cell

    def _is_site_overlapping(self, new_site, all_sites, inverses, inv):
        """Return True (and warn) if new_site overlaps with a site in all_sites."""
        diffs1 = np.abs(new_site - all_sites)
        diffs2 = np.abs(diffs1 - 1)
        mask = np.all(np.logical_or(diffs1 <= _Reader.equiv_site_tol,
                                    diffs2 <= _Reader.equiv_site_tol),
                        axis=-1)


        if np.any(mask):

            where_equal = np.argwhere(mask).flatten()
            for ind in where_equal:
                if inverses[ind] == inv:
                    pass
                else:
                    warnings.warn(
                        f'{self.current_name} has equivalent positions {inverses[ind]} and {inv}')
            return True
        else:
            return False

    def _validate_sites(self, asym_unit, asym_is_disordered):

        site_diffs1 = np.abs(asym_unit[:, None] - asym_unit)
        site_diffs2 = np.abs(site_diffs1 - 1)
        overlapping = np.triu(np.all(
            (site_diffs1 <= _Reader.equiv_site_tol) |
            (site_diffs2 <= _Reader.equiv_site_tol),
            axis=-1), 1)

        if self.disorder == 'all_sites':
            for i, j in np.argwhere(overlapping):
                if asym_is_disordered[i] or asym_is_disordered[j]:
                    overlapping[i, j] = False

        if overlapping.any():
            warnings.warn(
                f'{self.current_name} may have overlapping sites; duplicates will be removed')
            keep_sites = ~overlapping.any(0)
            return keep_sites

    def _has_no_valid_sites(self, motif):
        if motif.shape[0] == 0:
            warnings.warn(
                f'Skipping {self.current_name} as there are no sites with coordinates')
            return True
        return False

    def _construct_periodic_set(self, asym_unit, asym_symbols, sitesym, cell, **kwargs):
        """Asym motif + symbols + sitesym + cell (+kwargs) --> PeriodicSet"""
        frac_motif, asym_inds, multiplicities, inverses = self.expand(asym_unit, sitesym)
        full_types = [asym_symbols[i] for i in inverses]
        motif = frac_motif @ cell

        tags = {
            'name': self.current_name,
            'asymmetric_unit': asym_inds,
            'wyckoff_multiplicities': multiplicities,
            'types': full_types,
            **kwargs
        }

        if self.current_filename:
            tags['filename'] = self.current_filename

        return PeriodicSet(motif, cell, **tags)

    def expand(self, asym_unit: np.ndarray, sitesym: Sequence[str]) -> Tuple[np.ndarray, ...]:
        """
        Asymmetric unit's fractional coords + sitesyms (as strings)
        -->
        frac motif, asym unit inds, multiplicities, inverses
        """

        rotations, translations = ase.spacegroup.spacegroup.parse_sitesym(sitesym)
        all_sites = []
        asym_inds = [0]
        multiplicities = []
        inverses = []

        for inv, site in enumerate(asym_unit):
            multiplicity = 0

            for rot, trans in zip(rotations, translations):
                site_ = np.mod(np.dot(rot, site) + trans, 1)

                if not all_sites:
                    all_sites.append(site_)
                    inverses.append(inv)
                    multiplicity += 1
                    continue

                if not self._is_site_overlapping(site_, all_sites, inverses, inv):
                    all_sites.append(site_)
                    inverses.append(inv)
                    multiplicity += 1

            if multiplicity > 0:
                multiplicities.append(multiplicity)
                asym_inds.append(len(all_sites))

        frac_motif = np.array(all_sites)
        asym_inds = np.array(asym_inds[:-1])
        multiplicities = np.array(multiplicities)
        return frac_motif, asym_inds, multiplicities, inverses


def _atom_has_disorder(label, occupancy):
    return label.endswith('?') or (np.isscalar(occupancy) and occupancy < 1)

def _heaviest_component(molecule):
    """Heaviest component (removes all but the heaviest component of the asym unit).
    Intended for removing solvents. Probably doesn't play well with disorder"""
    component_weights = []
    for component in molecule.components:
        weight = 0
        for a in component.atoms:
            if isinstance(a.atomic_weight, (float, int)):
                if isinstance(a.occupancy, (float, int)):
                    weight += a.occupancy * a.atomic_weight
                else:
                    weight += a.atomic_weight
        component_weights.append(weight)
    largest_component_arg = np.argmax(np.array(component_weights))
    molecule = molecule.components[largest_component_arg]
    return molecule

def _validate_extract_data(extract_data):
    if not isinstance(extract_data, dict):
        raise ValueError('extract_data must be a dict of callables')
    for key in extract_data:
        if not callable(extract_data[key]):
            raise ValueError('extract_data must be a dict of callables')
        if key in _Reader.reserved_tags:
            raise ValueError(f'extract_data includes reserved key {key}')


1			"""Contains base reader class for the io module.
			0 ignored issues – show Coding Style introduced 2022-04-14 08:53 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
2			This class implements the converters from CifBlock, Entry to PeriodicSets.
3			"""
4
5			import warnings
6			from typing import Callable, Iterable, Sequence, Tuple
7
8			import numpy as np
9			import ase.spacegroup.spacegroup # parse_sitesym
10			import ase.io.cif
11
12			from .periodicset import PeriodicSet
13			from .utils import cellpar_to_cell
14
15
16			class _Reader:
			0 ignored issues – show best-practice introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report Too many instance attributes (10/7) Loading history...
17			"""Base Reader class. Contains parsers for converting ase CifBlock
18			and ccdc Entry objects to PeriodicSets.
19
20			Intended use:
21
22			First make a new method for _Reader converting object to PeriodicSet
23			(e.g. named _X_to_PSet). Then make this class outline:
24
25			class XReader(_Reader):
26			def __init__(self, ..., **kwargs):
27
28			super().__init__(**kwargs)
29
30			# setup and checks
31
32			# make 'iterable' which yields objects to be converted (e.g. CIFBlock, Entry)
33
34			# set self._generator like this
35			self._generator = self._read(iterable, self._X_to_PSet)
36			"""
37
38			disorder_options = {'skip', 'ordered_sites', 'all_sites'}
39			reserved_tags = {
40			'motif',
41			'cell',
42			'name',
43			'asymmetric_unit',
44			'wyckoff_multiplicities',
45			'types',
46			'filename',}
47			atom_site_fract_tags = [
48			'_atom_site_fract_x',
49			'_atom_site_fract_y',
50			'_atom_site_fract_z',]
51			atom_site_cartn_tags = [
52			'_atom_site_cartn_x',
53			'_atom_site_cartn_y',
54			'_atom_site_cartn_z',]
55			symop_tags = [
56			'_space_group_symop_operation_xyz',
57			'_space_group_symop.operation_xyz',
58			'_symmetry_equiv_pos_as_xyz',]
59
60			equiv_site_tol = 1e-3
61
62			def __init__(
			0 ignored issues – show best-practice introduced 2022-04-14 08:53 UTC by Report Bug Copy Issue Report Too many arguments (7/5) Loading history...
63			self,
64			remove_hydrogens=False,
65			disorder='skip',
66			heaviest_component=False,
67			show_warnings=True,
68			extract_data=None,
69			include_if=None):
70
71			if disorder not in _Reader.disorder_options:
72			raise ValueError(f'disorder parameter {disorder} must be one of {_Reader.disorder_options}')
			0 ignored issues – show Coding Style introduced 2022-04-14 13:33 UTC by Report Bug Copy Issue Report This line is too long as per the coding-style (104/100). This check looks for lines that are too long. You can specify the maximum line length. Loading history...
73
74			if extract_data is None:
75			self.extract_data = {}
76			else:
77			_validate_extract_data(extract_data)
78			self.extract_data = extract_data
79
80			if include_if is None:
81			self.include_if = ()
82			elif not all(callable(func) for func in include_if):
83			raise ValueError('include_if must be a list of callables')
84			else:
85			self.include_if = include_if
86
			0 ignored issues – show Coding Style introduced 2022-04-14 08:53 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
87			self.remove_hydrogens = remove_hydrogens
88			self.disorder = disorder
89			self.heaviest_component = heaviest_component
90			self.show_warnings = show_warnings
91			self.current_name = None
92			self.current_filename = None
93			self._generator = []
94
95			def __iter__(self):
96			yield from self._generator
97
98			def read_one(self):
99			"""Read the next (or first) item."""
100			return next(iter(self._generator))
101
102			def _map(self, func: Callable, iterable: Iterable) -> Iterable[PeriodicSet]:
103			"""Iterates over iterable, passing items through parser and yielding the
			0 ignored issues – show Coding Style introduced 2022-04-14 08:53 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
104			result if it is not None. Applies warning and include_if filter.
105			"""
106
107			with warnings.catch_warnings():
108			if not self.show_warnings:
109			warnings.simplefilter('ignore')
110			for item in iterable:
111			if all(check(item) for check in self.include_if):
112			res = func(item)
113			if res is not None:
114			yield res
115
116			def _cifblock_to_periodicset(self, block) -> PeriodicSet:
117			"""ase.io.cif.CIFBlock --> PeriodicSet. Returns None for a "bad" set."""
118
119			self.current_name = block.name
120			asym_unit, asym_symbols, sitesym, cell = self._cifblock_to_asym_unit(block)
121
122			# indices of sites to remove
123			remove = []
124			if self.remove_hydrogens:
125			remove.extend((i for i, sym in enumerate(asym_symbols) if sym in 'HD'))
			0 ignored issues – show Comprehensibility Best Practice introduced 2022-04-14 08:53 UTC by Report Bug Copy Issue Report The variable `i` does not seem to be defined. Loading history...
126
127			# find disordered sites
128			asym_is_disordered = []
129			occupancies = block.get('_atom_site_occupancy')
130			labels = block.get('_atom_site_label')
131			if occupancies is not None:
132			disordered = [] # indices where there is disorder
133			for i, (occ, label) in enumerate(zip(occupancies, labels)):
134			if _atom_has_disorder(label, occ):
135			if i not in remove:
136			disordered.append(i)
137			asym_is_disordered.append(True)
138			else:
139			asym_is_disordered.append(False)
140
141			if self.disorder == 'skip' and len(disordered) > 0:
142			warnings.warn(f'Skipping {self.current_name} as structure is disordered')
143			return None
144
145			if self.disorder == 'ordered_sites':
146			remove.extend(disordered)
147
148			# remove sites
149			asym_unit = np.delete(asym_unit, remove, axis=0)
150			asym_symbols = [s for i, s in enumerate(asym_symbols) if i not in remove]
151			asym_is_disordered = [v for i, v in enumerate(asym_is_disordered) if i not in remove]
152
153			keep_sites = self._validate_sites(asym_unit, asym_is_disordered)
154			if keep_sites is not None:
155			asym_unit = asym_unit[keep_sites]
156			asym_symbols = [sym for sym, keep in zip(asym_symbols, keep_sites) if keep]
157
158			if self._has_no_valid_sites(asym_unit):
159			return None
160
161			data = {key: func(block) for key, func in self.extract_data.items()}
162			periodic_set = self._construct_periodic_set(asym_unit, asym_symbols, sitesym, cell, **data)
163			return periodic_set
164
165			def _cifblock_to_asym_unit(self, block):
166			"""ase.io.cif.CIFBlock -->
			0 ignored issues – show Coding Style introduced 2022-04-14 13:51 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
167			asymmetric unit (frac coords), asym_symbols, cell, symops (as strings)"""
168
			0 ignored issues – show Coding Style introduced 2022-04-15 15:50 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
169			cell = block.get_cell().array
170			asym_unit = [block.get(name) for name in _Reader.atom_site_fract_tags]
171			if None in asym_unit:
172			asym_motif = [block.get(name) for name in _Reader.atom_site_cartn_tags]
173			if None in asym_motif:
174			warnings.warn(f'Skipping {self.current_name} as coordinates were not found')
175			return None
176			asym_unit = np.array(asym_motif) @ np.linalg.inv(cell)
177			asym_unit = np.mod(np.array(asym_unit).T, 1)
178
			0 ignored issues – show Coding Style introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
179			try:
180			asym_symbols = block.get_symbols()
181			except ase.io.cif.NoStructureData as _:
182			asym_symbols = ['Unknown' for _ in range(len(asym_unit))]
183
			0 ignored issues – show Coding Style introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
184			sitesym = ['x,y,z', ]
185			for tag in _Reader.symop_tags:
186			if tag in block:
187			sitesym = block[tag]
188			break
189
190			if isinstance(sitesym, str):
191			sitesym = [sitesym]
192
			0 ignored issues – show Coding Style introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
193			return asym_unit, asym_symbols, sitesym, cell
194
195			def _entry_to_periodicset(self, entry) -> PeriodicSet:
196			"""ccdc.entry.Entry --> PeriodicSet. Returns None for a "bad" set."""
197
198			self.current_name = entry.identifier
199			crystal = entry.crystal
200
201			if not entry.has_3d_structure:
202			warnings.warn(f'Skipping {self.current_name} as entry has no 3D structure')
203			return None
204
205			# first disorder check, if skipping. If occ == 1 for all atoms but the entry
206			# or crystal is listed as having disorder, skip (can't know where disorder is).
207			# If occ != 1 for any atoms, we wait to see if we remove them before skipping.
208			molecule = crystal.disordered_molecule
209			if self.disorder == 'ordered_sites':
210			molecule.remove_atoms(a for a in molecule.atoms if a.label.endswith('?'))
211
212			may_have_disorder = False
213			if self.disorder == 'skip':
214			for a in molecule.atoms:
215			occ = a.occupancy
216			if _atom_has_disorder(a.label, occ):
217			may_have_disorder = True
218			break
219
220			if not may_have_disorder:
221			if crystal.has_disorder or entry.has_disorder:
222			warnings.warn(f'Skipping {self.current_name} as structure is disordered')
223			return None
224
225			# make same as cifblock version??
226			if self.remove_hydrogens:
227			molecule.remove_atoms(a for a in molecule.atoms if a.atomic_symbol in 'HD')
228
229			if self.heaviest_component and len(molecule.components) > 1:
230			molecule = _heaviest_component(molecule)
231
232			crystal.molecule = molecule
233
234			# by here all atoms to be removed have been (except via ordered_sites).
235			# If disorder == 'skip' and there were atom(s) with occ < 1 found
236			# eariler, we check if all such atoms were removed. If not, skip.
237			if self.disorder == 'skip' and may_have_disorder:
238			for a in crystal.disordered_molecule.atoms:
239			occ = a.occupancy
240			if _atom_has_disorder(a.label, occ):
241			warnings.warn(f'Skipping {self.current_name} as structure is disordered')
242			return None
243
244			# if disorder is all_sites, we need to know where disorder is to ignore overlaps
245			asym_is_disordered = [] # True/False list same length as asym unit
246			if self.disorder == 'all_sites':
247			for a in crystal.asymmetric_unit_molecule.atoms:
248			occ = a.occupancy
249			if _atom_has_disorder(a.label, occ):
250			asym_is_disordered.append(True)
251			else:
252			asym_is_disordered.append(False)
253
254			# check all atoms have coords. option/default remove unknown sites?
255			if not molecule.all_atoms_have_sites or \
256			any(a.fractional_coordinates is None for a in molecule.atoms):
257			warnings.warn(f'Skipping {self.current_name} as some atoms do not have sites')
258			return None
259
260			asym_unit, asym_symbols, sitesym, cell = self._crystal_to_asym_unit(crystal)
261
262			# remove overlapping sites, check sites exist
263			keep_sites = self._validate_sites(asym_unit, asym_is_disordered)
264			if keep_sites is not None:
265			asym_unit = asym_unit[keep_sites]
266			asym_symbols = [sym for sym, keep in zip(asym_symbols, keep_sites) if keep]
267
268			if self._has_no_valid_sites(asym_unit):
269			return None
270
271			entry.crystal.molecule = crystal.disordered_molecule
272			data = {key: func(entry) for key, func in self.extract_data.items()}
273			periodic_set = self._construct_periodic_set(asym_unit, asym_symbols, sitesym, cell, **data)
274			return periodic_set
275
276			def _crystal_to_asym_unit(self, crystal):
			0 ignored issues – show Coding Style introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report This method could be written as a function/class method. If a method does not access any attributes of the class, it could also be implemented as a function or static method. This can help improve readability. For example class Foo: def some_method(self, x, y): return x + y; could be written as class Foo: @classmethod def some_method(cls, x, y): return x + y; Loading history...
277			"""ase.io.cif.CIFBlock -->
278			asymmetric unit (frac coords), asym_symbols, symops, cell"""
279
			0 ignored issues – show Coding Style introduced 2022-04-15 17:02 UTC by Report Bug Copy Issue Report Trailing whitespace Loading history...
280			asym_atoms = crystal.asymmetric_unit_molecule.atoms
281			asym_unit = np.array([tuple(a.fractional_coordinates) for a in asym_atoms])
282			asym_unit = np.mod(asym_unit, 1)
283			asym_symbols = [a.atomic_symbol for a in asym_atoms]
284			cell = cellpar_to_cell(crystal.cell_lengths, crystal.cell_angles)
285			sitesym = crystal.symmetry_operators
286			if not sitesym:
287			sitesym = ['x,y,z', ]
288			return asym_unit, asym_symbols, sitesym, cell
289
290			def _is_site_overlapping(self, new_site, all_sites, inverses, inv):
291			"""Return True (and warn) if new_site overlaps with a site in all_sites."""
292			diffs1 = np.abs(new_site - all_sites)
293			diffs2 = np.abs(diffs1 - 1)
294			mask = np.all(np.logical_or(diffs1 <= _Reader.equiv_site_tol,
295			diffs2 <= _Reader.equiv_site_tol),
296			axis=-1)
			0 ignored issues – show Coding Style introduced 2022-04-14 12:12 UTC by Report Bug Copy Issue Report Wrong continued indentation (remove 2 spaces). Loading history...
297
298			if np.any(mask):
			0 ignored issues – show unused-code introduced 2022-04-14 12:12 UTC by Report Bug Copy Issue Report Unnecessary "else" after "return" Loading history...
299			where_equal = np.argwhere(mask).flatten()
300			for ind in where_equal:
301			if inverses[ind] == inv:
302			pass
303			else:
304			warnings.warn(
305			f'{self.current_name} has equivalent positions {inverses[ind]} and {inv}')
306			return True
307			else:
308			return False
309
310			def _validate_sites(self, asym_unit, asym_is_disordered):
			0 ignored issues – show Unused Code introduced 2022-04-14 13:03 UTC by Report Bug Copy Issue Report Either all return statements in a function should return an expression, or none of them should. Loading history...
311			site_diffs1 = np.abs(asym_unit[:, None] - asym_unit)
312			site_diffs2 = np.abs(site_diffs1 - 1)
313			overlapping = np.triu(np.all(
314			(site_diffs1 <= _Reader.equiv_site_tol) \|
315			(site_diffs2 <= _Reader.equiv_site_tol),
316			axis=-1), 1)
317
318			if self.disorder == 'all_sites':
319			for i, j in np.argwhere(overlapping):
320			if asym_is_disordered[i] or asym_is_disordered[j]:
321			overlapping[i, j] = False
322
323			if overlapping.any():
324			warnings.warn(
325			f'{self.current_name} may have overlapping sites; duplicates will be removed')
326			keep_sites = ~overlapping.any(0)
327			return keep_sites
328
329			def _has_no_valid_sites(self, motif):
330			if motif.shape[0] == 0:
331			warnings.warn(
332			f'Skipping {self.current_name} as there are no sites with coordinates')
333			return True
334			return False
335
336			def _construct_periodic_set(self, asym_unit, asym_symbols, sitesym, cell, **kwargs):
337			"""Asym motif + symbols + sitesym + cell (+kwargs) --> PeriodicSet"""
338			frac_motif, asym_inds, multiplicities, inverses = self.expand(asym_unit, sitesym)
339			full_types = [asym_symbols[i] for i in inverses]
340			motif = frac_motif @ cell
341
342			tags = {
343			'name': self.current_name,
344			'asymmetric_unit': asym_inds,
345			'wyckoff_multiplicities': multiplicities,
346			'types': full_types,
347			**kwargs
348			}
349
350			if self.current_filename:
351			tags['filename'] = self.current_filename
352
353			return PeriodicSet(motif, cell, **tags)
354
355			def expand(self, asym_unit: np.ndarray, sitesym: Sequence[str]) -> Tuple[np.ndarray, ...]:
356			"""
357			Asymmetric unit's fractional coords + sitesyms (as strings)
358			-->
359			frac motif, asym unit inds, multiplicities, inverses
360			"""
361
362			rotations, translations = ase.spacegroup.spacegroup.parse_sitesym(sitesym)
363			all_sites = []
364			asym_inds = [0]
365			multiplicities = []
366			inverses = []
367
368			for inv, site in enumerate(asym_unit):
369			multiplicity = 0
370
371			for rot, trans in zip(rotations, translations):
372			site_ = np.mod(np.dot(rot, site) + trans, 1)
373
374			if not all_sites:
375			all_sites.append(site_)
376			inverses.append(inv)
377			multiplicity += 1
378			continue
379
380			if not self._is_site_overlapping(site_, all_sites, inverses, inv):
381			all_sites.append(site_)
382			inverses.append(inv)
383			multiplicity += 1
384
385			if multiplicity > 0:
386			multiplicities.append(multiplicity)
387			asym_inds.append(len(all_sites))
388
389			frac_motif = np.array(all_sites)
390			asym_inds = np.array(asym_inds[:-1])
391			multiplicities = np.array(multiplicities)
392			return frac_motif, asym_inds, multiplicities, inverses
393
394
395			def _atom_has_disorder(label, occupancy):
396			return label.endswith('?') or (np.isscalar(occupancy) and occupancy < 1)
397
398			def _heaviest_component(molecule):
399			"""Heaviest component (removes all but the heaviest component of the asym unit).
400			Intended for removing solvents. Probably doesn't play well with disorder"""
401			component_weights = []
402			for component in molecule.components:
403			weight = 0
404			for a in component.atoms:
405			if isinstance(a.atomic_weight, (float, int)):
406			if isinstance(a.occupancy, (float, int)):
407			weight += a.occupancy * a.atomic_weight
408			else:
409			weight += a.atomic_weight
410			component_weights.append(weight)
411			largest_component_arg = np.argmax(np.array(component_weights))
412			molecule = molecule.components[largest_component_arg]
413			return molecule
414
415			def _validate_extract_data(extract_data):
416			if not isinstance(extract_data, dict):
417			raise ValueError('extract_data must be a dict of callables')
418			for key in extract_data:
419			if not callable(extract_data[key]):
420			raise ValueError('extract_data must be a dict of callables')
421			if key in _Reader.reserved_tags:
422			raise ValueError(f'extract_data includes reserved key {key}')
423

dwiddo / average-minimum-distance

Push — master ( ca5a49...2d0e3f )

amd._reader._Reader._construct_periodic_set() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like