From 14d20b1164d4ba022e1772702e4ae4cec60912cf Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Tue, 9 Jul 2024 18:33:41 +0200 Subject: [PATCH 01/16] Adjusted dependency path in pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index da497aba..44b903cc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,6 +31,7 @@ license = { file = "LICENSE" } dependencies = [ "nomad-lab>=1.3.0", "matid>=2.0.0.dev2", + "nomad-simulations@file:///home/bmohr/software/nomad-simulations", ] [project.urls] @@ -136,5 +137,4 @@ where = ["src"] write_to = "src/nomad_simulations/_version.py" [project.entry-points.'nomad.plugin'] - nomad_simulations_plugin = "nomad_simulations.schema_packages:nomad_simulations_plugin" From b30c9786ca0b1e4b5c36e6f09613206c29193afe Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Tue, 9 Jul 2024 18:37:57 +0200 Subject: [PATCH 02/16] reverted pyproject.toml --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 44b903cc..ff3ceb38 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,6 @@ license = { file = "LICENSE" } dependencies = [ "nomad-lab>=1.3.0", "matid>=2.0.0.dev2", - "nomad-simulations@file:///home/bmohr/software/nomad-simulations", ] [project.urls] @@ -137,4 +136,5 @@ where = ["src"] write_to = "src/nomad_simulations/_version.py" [project.entry-points.'nomad.plugin'] -nomad_simulations_plugin = "nomad_simulations.schema_packages:nomad_simulations_plugin" + +nomad_simulations_plugin = "nomad_simulations.schema_packages:nomad_simulations_plugin" \ No newline at end of file From 932add36286e83ee77cad599a54da4a6ecac84ec Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Thu, 11 Jul 2024 15:05:47 +0200 Subject: [PATCH 03/16] test --- pyproject.toml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index ff3ceb38..0cd48773 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,11 @@ maintainers = [ { name = "Joseph F. Rudzinski", email = "joseph.rudzinski@physik.hu-berlin.de" } ] license = { file = "LICENSE" } +# dependencies = [ +# "nomad-lab>=1.3.0", +# "matid>=2.0.0.dev2", +# "nomad-simulations@file:///home/bmohr/software/nomad-simulations", +# ] dependencies = [ "nomad-lab>=1.3.0", "matid>=2.0.0.dev2", @@ -137,4 +142,4 @@ write_to = "src/nomad_simulations/_version.py" [project.entry-points.'nomad.plugin'] -nomad_simulations_plugin = "nomad_simulations.schema_packages:nomad_simulations_plugin" \ No newline at end of file +nomad_simulations_plugin = "nomad_simulations.schema_packages:nomad_simulations_plugin" From 4e2b5596f060405004f8e729d707c4e9ac1dc790 Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Tue, 8 Oct 2024 11:03:01 +0200 Subject: [PATCH 04/16] Adding support for coarse-grained models --- .../schema_packages/particles_state.py | 2132 +++++++++++++++++ 1 file changed, 2132 insertions(+) create mode 100644 src/nomad_simulations/schema_packages/particles_state.py diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py new file mode 100644 index 00000000..458a00d3 --- /dev/null +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -0,0 +1,2132 @@ +from typing import TYPE_CHECKING, Any, Optional, Union + +import numpy as np +import pint +from nomad.datamodel.data import ArchiveSection +from nomad.datamodel.metainfo.annotations import ELNAnnotation +from nomad.datamodel.metainfo.basesections import Entity +from nomad.metainfo import MEnum, Quantity, SubSection +from nomad.units import ureg + +if TYPE_CHECKING: + from nomad.datamodel.datamodel import EntryArchive + from nomad.metainfo import Context, Section + from structlog.stdlib import BoundLogger + + +class Particle: + """Particle object. + + The Particle object can represent an isolated compound, or a + periodically repeated structure. It has a unit cell and + there may be periodic boundary conditions along any of the three + unit cell axes. + Information about the particles (types and position) is + stored in ndarrays. Optionally, there can be information about + type_shapes, mass, charge, diameter, body, moment_inertia, orientation, + angular momenta and image. + + In order to calculate energies, forces and stresses, a calculator + object has to attached to the atoms object. + + Parameters: + + types: str (formula) or list of str + Can be a string formula, a list of symbols or a list of + Atom objects. Examples: 'H2O', 'COPt12', ['H', 'H', 'O'], + [Atom('Ne', (x, y, z)), ...]. + positions: list of xyz-positions + Atomic positions. Anything that can be converted to an + ndarray of shape (n, 3) will do: [(x1,y1,z1), (x2,y2,z2), + ...]. + scaled_positions: list of scaled-positions + Like positions, but given in units of the unit cell. + Can not be set at the same time as positions. + numbers: list of int + Atomic numbers (use only one of symbols/numbers). + tags: list of int + Special purpose tags. + momenta: list of xyz-momenta + Momenta for all atoms. + masses: list of float + Atomic masses in atomic units. + magmoms: list of float or list of xyz-values + Magnetic moments. Can be either a single value for each atom + for collinear calculations or three numbers for each atom for + non-collinear calculations. + charges: list of float + Initial atomic charges. + cell: 3x3 matrix or length 3 or 6 vector + Unit cell vectors. Can also be given as just three + numbers for orthorhombic cells, or 6 numbers, where + first three are lengths of unit cell vectors, and the + other three are angles between them (in degrees), in following order: + [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]. + First vector will lie in x-direction, second in xy-plane, + and the third one in z-positive subspace. + Default value: [0, 0, 0]. + celldisp: Vector + Unit cell displacement vector. To visualize a displaced cell + around the center of mass of a Systems of atoms. Default value + = (0,0,0) + pbc: one or three bool + Periodic boundary conditions flags. Examples: True, + False, 0, 1, (1, 1, 0), (True, False, False). Default + value: False. + constraint: constraint object(s) + Used for applying one or more constraints during structure + optimization. + calculator: calculator object + Used to attach a calculator for calculating energies and atomic + forces. + info: dict of key-value pairs + Dictionary of key-value pairs with additional information + about the system. The following keys may be used by ase: + + - spacegroup: Spacegroup instance + - unit_cell: 'conventional' | 'primitive' | int | 3 ints + - adsorbate_info: Information about special adsorption sites + + Items in the info attribute survives copy and slicing and can + be stored in and retrieved from trajectory files given that the + key is a string, the value is JSON-compatible and, if the value is a + user-defined object, its base class is importable. One should + not make any assumptions about the existence of keys. + + Examples: + + These three are equivalent: + + >>> d = 1.104 # N2 bondlength + >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)]) + >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)]) + >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))]) + + FCC gold: + + >>> a = 4.05 # Gold lattice constant + >>> b = a / 2 + >>> fcc = Atoms('Au', + ... cell=[(0, b, b), (b, 0, b), (b, b, 0)], + ... pbc=True) + + Hydrogen wire: + + >>> d = 0.9 # H-H distance + >>> h = Atoms('H', positions=[(0, 0, 0)], + ... cell=(d, 0, 0), + ... pbc=(1, 0, 0)) + """ + + ase_objtype = 'atoms' # For JSONability + + def __init__( + self, + symbols=None, + positions=None, + numbers=None, + tags=None, + momenta=None, + masses=None, + magmoms=None, + charges=None, + scaled_positions=None, + cell=None, + pbc=None, + celldisp=None, + constraint=None, + calculator=None, + info=None, + velocities=None, + ): + self._cellobj = Cell.new() + self._pbc = np.zeros(3, bool) + + atoms = None + + if hasattr(symbols, 'get_positions'): + atoms = symbols + symbols = None + elif ( + isinstance(symbols, (list, tuple)) + and len(symbols) > 0 + and isinstance(symbols[0], Atom) + ): + # Get data from a list or tuple of Atom objects: + data = [ + [atom.get_raw(name) for atom in symbols] + for name in [ + 'position', + 'number', + 'tag', + 'momentum', + 'mass', + 'magmom', + 'charge', + ] + ] + atoms = self.__class__(None, *data) + symbols = None + + if atoms is not None: + # Get data from another Atoms object: + if scaled_positions is not None: + raise NotImplementedError + if symbols is None and numbers is None: + numbers = atoms.get_atomic_numbers() + if positions is None: + positions = atoms.get_positions() + if tags is None and atoms.has('tags'): + tags = atoms.get_tags() + if momenta is None and atoms.has('momenta'): + momenta = atoms.get_momenta() + if magmoms is None and atoms.has('initial_magmoms'): + magmoms = atoms.get_initial_magnetic_moments() + if masses is None and atoms.has('masses'): + masses = atoms.get_masses() + if charges is None and atoms.has('initial_charges'): + charges = atoms.get_initial_charges() + if cell is None: + cell = atoms.get_cell() + if celldisp is None: + celldisp = atoms.get_celldisp() + if pbc is None: + pbc = atoms.get_pbc() + if constraint is None: + constraint = [c.copy() for c in atoms.constraints] + if calculator is None: + calculator = atoms.calc + if info is None: + info = copy.deepcopy(atoms.info) + + self.arrays = {} + + if symbols is None: + if numbers is None: + if positions is not None: + natoms = len(positions) + elif scaled_positions is not None: + natoms = len(scaled_positions) + else: + natoms = 0 + numbers = np.zeros(natoms, int) + self.new_array('numbers', numbers, int) + else: + if numbers is not None: + raise TypeError('Use only one of "symbols" and "numbers".') + else: + self.new_array('numbers', symbols2numbers(symbols), int) + + if self.numbers.ndim != 1: + raise ValueError('"numbers" must be 1-dimensional.') + + if cell is None: + cell = np.zeros((3, 3)) + self.set_cell(cell) + + if celldisp is None: + celldisp = np.zeros(shape=(3, 1)) + self.set_celldisp(celldisp) + + if positions is None: + if scaled_positions is None: + positions = np.zeros((len(self.arrays['numbers']), 3)) + else: + assert self.cell.rank == 3 + positions = np.dot(scaled_positions, self.cell) + else: + if scaled_positions is not None: + raise TypeError('Use only one of "symbols" and "numbers".') + self.new_array('positions', positions, float, (3,)) + + self.set_constraint(constraint) + self.set_tags(default(tags, 0)) + self.set_masses(default(masses, None)) + self.set_initial_magnetic_moments(default(magmoms, 0.0)) + self.set_initial_charges(default(charges, 0.0)) + if pbc is None: + pbc = False + self.set_pbc(pbc) + self.set_momenta(default(momenta, (0.0, 0.0, 0.0)), apply_constraint=False) + + if velocities is not None: + if momenta is None: + self.set_velocities(velocities) + else: + raise TypeError('Use only one of "momenta" and "velocities".') + + if info is None: + self.info = {} + else: + self.info = dict(info) + + self.calc = calculator + + @property + def symbols(self): + """Get chemical symbols as a :class:`ase.symbols.Symbols` object. + + The object works like ``atoms.numbers`` except its values + are strings. It supports in-place editing.""" + return Symbols(self.numbers) + + @symbols.setter + def symbols(self, obj): + new_symbols = Symbols.fromsymbols(obj) + self.numbers[:] = new_symbols.numbers + + @deprecated(DeprecationWarning('Please use atoms.calc = calc')) + def set_calculator(self, calc=None): + """Attach calculator object. + + Please use the equivalent atoms.calc = calc instead of this + method.""" + self.calc = calc + + @deprecated(DeprecationWarning('Please use atoms.calc')) + def get_calculator(self): + """Get currently attached calculator object. + + Please use the equivalent atoms.calc instead of + atoms.get_calculator().""" + return self.calc + + @property + def calc(self): + """Calculator object.""" + return self._calc + + @calc.setter + def calc(self, calc): + self._calc = calc + if hasattr(calc, 'set_atoms'): + calc.set_atoms(self) + + @calc.deleter # type: ignore + @deprecated(DeprecationWarning('Please use atoms.calc = None')) + def calc(self): + self._calc = None + + @property # type: ignore + @deprecated('Please use atoms.cell.rank instead') + def number_of_lattice_vectors(self): + """Number of (non-zero) lattice vectors.""" + return self.cell.rank + + def set_constraint(self, constraint=None): + """Apply one or more constrains. + + The *constraint* argument must be one constraint object or a + list of constraint objects.""" + if constraint is None: + self._constraints = [] + else: + if isinstance(constraint, list): + self._constraints = constraint + elif isinstance(constraint, tuple): + self._constraints = list(constraint) + else: + self._constraints = [constraint] + + def _get_constraints(self): + return self._constraints + + def _del_constraints(self): + self._constraints = [] + + constraints = property( + _get_constraints, set_constraint, _del_constraints, 'Constraints of the atoms.' + ) + + def set_cell(self, cell, scale_atoms=False, apply_constraint=True): + """Set unit cell vectors. + + Parameters: + + cell: 3x3 matrix or length 3 or 6 vector + Unit cell. A 3x3 matrix (the three unit cell vectors) or + just three numbers for an orthorhombic cell. Another option is + 6 numbers, which describes unit cell with lengths of unit cell + vectors and with angles between them (in degrees), in following + order: [len(a), len(b), len(c), angle(b,c), angle(a,c), + angle(a,b)]. First vector will lie in x-direction, second in + xy-plane, and the third one in z-positive subspace. + scale_atoms: bool + Fix atomic positions or move atoms with the unit cell? + Default behavior is to *not* move the atoms (scale_atoms=False). + apply_constraint: bool + Whether to apply constraints to the given cell. + + Examples: + + Two equivalent ways to define an orthorhombic cell: + + >>> atoms = Atoms('He') + >>> a, b, c = 7, 7.5, 8 + >>> atoms.set_cell([a, b, c]) + >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)]) + + FCC unit cell: + + >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)]) + + Hexagonal unit cell: + + >>> atoms.set_cell([a, a, c, 90, 90, 120]) + + Rhombohedral unit cell: + + >>> alpha = 77 + >>> atoms.set_cell([a, a, a, alpha, alpha, alpha]) + """ + + # Override pbcs if and only if given a Cell object: + cell = Cell.new(cell) + + # XXX not working well during initialize due to missing _constraints + if apply_constraint and hasattr(self, '_constraints'): + for constraint in self.constraints: + if hasattr(constraint, 'adjust_cell'): + constraint.adjust_cell(self, cell) + + if scale_atoms: + M = np.linalg.solve(self.cell.complete(), cell.complete()) + self.positions[:] = np.dot(self.positions, M) + + self.cell[:] = cell + + def set_celldisp(self, celldisp): + """Set the unit cell displacement vectors.""" + celldisp = np.array(celldisp, float) + self._celldisp = celldisp + + def get_celldisp(self): + """Get the unit cell displacement vectors.""" + return self._celldisp.copy() + + def get_cell(self, complete=False): + """Get the three unit cell vectors as a `class`:ase.cell.Cell` object. + + The Cell object resembles a 3x3 ndarray, and cell[i, j] + is the jth Cartesian coordinate of the ith cell vector.""" + if complete: + cell = self.cell.complete() + else: + cell = self.cell.copy() + + return cell + + @deprecated('Please use atoms.cell.cellpar() instead') + def get_cell_lengths_and_angles(self): + """Get unit cell parameters. Sequence of 6 numbers. + + First three are unit cell vector lengths and second three + are angles between them:: + + [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)] + + in degrees. + """ + return self.cell.cellpar() + + @deprecated('Please use atoms.cell.reciprocal()') + def get_reciprocal_cell(self): + """Get the three reciprocal lattice vectors as a 3x3 ndarray. + + Note that the commonly used factor of 2 pi for Fourier + transforms is not included here.""" + + return self.cell.reciprocal() + + @property + def pbc(self): + """Reference to pbc-flags for in-place manipulations.""" + return self._pbc + + @pbc.setter + def pbc(self, pbc): + self._pbc[:] = pbc + + def set_pbc(self, pbc): + """Set periodic boundary condition flags.""" + self.pbc = pbc + + def get_pbc(self): + """Get periodic boundary condition flags.""" + return self.pbc.copy() + + def new_array(self, name, a, dtype=None, shape=None): + """Add new array. + + If *shape* is not *None*, the shape of *a* will be checked.""" + + if dtype is not None: + a = np.array(a, dtype, order='C') + if len(a) == 0 and shape is not None: + a.shape = (-1,) + shape + else: + if not a.flags['C_CONTIGUOUS']: + a = np.ascontiguousarray(a) + else: + a = a.copy() + + if name in self.arrays: + raise RuntimeError(f'Array {name} already present') + + for b in self.arrays.values(): + if len(a) != len(b): + raise ValueError( + 'Array "%s" has wrong length: %d != %d.' % (name, len(a), len(b)) + ) + break + + if shape is not None and a.shape[1:] != shape: + raise ValueError( + 'Array "%s" has wrong shape %s != %s.' + % (name, a.shape, (a.shape[0:1] + shape)) + ) + + self.arrays[name] = a + + def get_array(self, name, copy=True): + """Get an array. + + Returns a copy unless the optional argument copy is false. + """ + if copy: + return self.arrays[name].copy() + else: + return self.arrays[name] + + def set_array(self, name, a, dtype=None, shape=None): + """Update array. + + If *shape* is not *None*, the shape of *a* will be checked. + If *a* is *None*, then the array is deleted.""" + + b = self.arrays.get(name) + if b is None: + if a is not None: + self.new_array(name, a, dtype, shape) + else: + if a is None: + del self.arrays[name] + else: + a = np.asarray(a) + if a.shape != b.shape: + raise ValueError( + 'Array "%s" has wrong shape %s != %s.' + % (name, a.shape, b.shape) + ) + b[:] = a + + def has(self, name): + """Check for existence of array. + + name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms', + 'initial_charges'.""" + # XXX extend has to calculator properties + return name in self.arrays + + def set_atomic_numbers(self, numbers): + """Set atomic numbers.""" + self.set_array('numbers', numbers, int, ()) + + def get_atomic_numbers(self): + """Get integer array of atomic numbers.""" + return self.arrays['numbers'].copy() + + def get_chemical_symbols(self): + """Get list of chemical symbol strings. + + Equivalent to ``list(atoms.symbols)``.""" + return list(self.symbols) + + def set_chemical_symbols(self, symbols): + """Set chemical symbols.""" + self.set_array('numbers', symbols2numbers(symbols), int, ()) + + def get_chemical_formula(self, mode='hill', empirical=False): + """Get the chemical formula as a string based on the chemical symbols. + + Parameters: + + mode: str + There are four different modes available: + + 'all': The list of chemical symbols are contracted to a string, + e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'. + + 'reduce': The same as 'all' where repeated elements are contracted + to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to + 'CH3OCH3'. + + 'hill': The list of chemical symbols are contracted to a string + following the Hill notation (alphabetical order with C and H + first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to + 'H2O4S'. This is default. + + 'metal': The list of chemical symbols (alphabetical metals, + and alphabetical non-metals) + + empirical, bool (optional, default=False) + Divide the symbol counts by their greatest common divisor to yield + an empirical formula. Only for mode `metal` and `hill`. + """ + return self.symbols.get_chemical_formula(mode, empirical) + + def set_tags(self, tags): + """Set tags for all atoms. If only one tag is supplied, it is + applied to all atoms.""" + if isinstance(tags, int): + tags = [tags] * len(self) + self.set_array('tags', tags, int, ()) + + def get_tags(self): + """Get integer array of tags.""" + if 'tags' in self.arrays: + return self.arrays['tags'].copy() + else: + return np.zeros(len(self), int) + + def set_momenta(self, momenta, apply_constraint=True): + """Set momenta.""" + if apply_constraint and len(self.constraints) > 0 and momenta is not None: + momenta = np.array(momenta) # modify a copy + for constraint in self.constraints: + if hasattr(constraint, 'adjust_momenta'): + constraint.adjust_momenta(self, momenta) + self.set_array('momenta', momenta, float, (3,)) + + def set_velocities(self, velocities): + """Set the momenta by specifying the velocities.""" + self.set_momenta(self.get_masses()[:, np.newaxis] * velocities) + + def get_momenta(self): + """Get array of momenta.""" + if 'momenta' in self.arrays: + return self.arrays['momenta'].copy() + else: + return np.zeros((len(self), 3)) + + def set_masses(self, masses='defaults'): + """Set atomic masses in atomic mass units. + + The array masses should contain a list of masses. In case + the masses argument is not given or for those elements of the + masses list that are None, standard values are set.""" + + if isinstance(masses, str): + if masses == 'defaults': + masses = atomic_masses[self.arrays['numbers']] + elif masses == 'most_common': + masses = atomic_masses_common[self.arrays['numbers']] + elif masses is None: + pass + elif not isinstance(masses, np.ndarray): + masses = list(masses) + for i, mass in enumerate(masses): + if mass is None: + masses[i] = atomic_masses[self.numbers[i]] + self.set_array('masses', masses, float, ()) + + def get_masses(self): + """Get array of masses in atomic mass units.""" + if 'masses' in self.arrays: + return self.arrays['masses'].copy() + else: + return atomic_masses[self.arrays['numbers']] + + def set_initial_magnetic_moments(self, magmoms=None): + """Set the initial magnetic moments. + + Use either one or three numbers for every atom (collinear + or non-collinear spins).""" + + if magmoms is None: + self.set_array('initial_magmoms', None) + else: + magmoms = np.asarray(magmoms) + self.set_array('initial_magmoms', magmoms, float, magmoms.shape[1:]) + + def get_initial_magnetic_moments(self): + """Get array of initial magnetic moments.""" + if 'initial_magmoms' in self.arrays: + return self.arrays['initial_magmoms'].copy() + else: + return np.zeros(len(self)) + + def get_magnetic_moments(self): + """Get calculated local magnetic moments.""" + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + return self._calc.get_magnetic_moments(self) + + def get_magnetic_moment(self): + """Get calculated total magnetic moment.""" + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + return self._calc.get_magnetic_moment(self) + + def set_initial_charges(self, charges=None): + """Set the initial charges.""" + + if charges is None: + self.set_array('initial_charges', None) + else: + self.set_array('initial_charges', charges, float, ()) + + def get_initial_charges(self): + """Get array of initial charges.""" + if 'initial_charges' in self.arrays: + return self.arrays['initial_charges'].copy() + else: + return np.zeros(len(self)) + + def get_charges(self): + """Get calculated charges.""" + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + try: + return self._calc.get_charges(self) + except AttributeError: + from ase.calculators.calculator import PropertyNotImplementedError + + raise PropertyNotImplementedError + + def set_positions(self, newpositions, apply_constraint=True): + """Set positions, honoring any constraints. To ignore constraints, + use *apply_constraint=False*.""" + if self.constraints and apply_constraint: + newpositions = np.array(newpositions, float) + for constraint in self.constraints: + constraint.adjust_positions(self, newpositions) + + self.set_array('positions', newpositions, shape=(3,)) + + def get_positions(self, wrap=False, **wrap_kw): + """Get array of positions. + + Parameters: + + wrap: bool + wrap atoms back to the cell before returning positions + wrap_kw: (keyword=value) pairs + optional keywords `pbc`, `center`, `pretty_translation`, `eps`, + see :func:`ase.geometry.wrap_positions` + """ + if wrap: + if 'pbc' not in wrap_kw: + wrap_kw['pbc'] = self.pbc + return wrap_positions(self.positions, self.cell, **wrap_kw) + else: + return self.arrays['positions'].copy() + + def get_potential_energy(self, force_consistent=False, apply_constraint=True): + """Calculate potential energy. + + Ask the attached calculator to calculate the potential energy and + apply constraints. Use *apply_constraint=False* to get the raw + forces. + + When supported by the calculator, either the energy extrapolated + to zero Kelvin or the energy consistent with the forces (the free + energy) can be returned. + """ + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + if force_consistent: + energy = self._calc.get_potential_energy( + self, force_consistent=force_consistent + ) + else: + energy = self._calc.get_potential_energy(self) + if apply_constraint: + for constraint in self.constraints: + if hasattr(constraint, 'adjust_potential_energy'): + energy += constraint.adjust_potential_energy(self) + return energy + + def get_properties(self, properties): + """This method is experimental; currently for internal use.""" + # XXX Something about constraints. + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + return self._calc.calculate_properties(self, properties) + + def get_potential_energies(self): + """Calculate the potential energies of all the atoms. + + Only available with calculators supporting per-atom energies + (e.g. classical potentials). + """ + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + return self._calc.get_potential_energies(self) + + def get_kinetic_energy(self): + """Get the kinetic energy.""" + momenta = self.arrays.get('momenta') + if momenta is None: + return 0.0 + return 0.5 * np.vdot(momenta, self.get_velocities()) + + def get_velocities(self): + """Get array of velocities.""" + momenta = self.get_momenta() + masses = self.get_masses() + return momenta / masses[:, np.newaxis] + + def get_total_energy(self): + """Get the total energy - potential plus kinetic energy.""" + return self.get_potential_energy() + self.get_kinetic_energy() + + def get_forces(self, apply_constraint=True, md=False): + """Calculate atomic forces. + + Ask the attached calculator to calculate the forces and apply + constraints. Use *apply_constraint=False* to get the raw + forces. + + For molecular dynamics (md=True) we don't apply the constraint + to the forces but to the momenta. When holonomic constraints for + rigid linear triatomic molecules are present, ask the constraints + to redistribute the forces within each triple defined in the + constraints (required for molecular dynamics with this type of + constraints).""" + + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + forces = self._calc.get_forces(self) + + if apply_constraint: + # We need a special md flag here because for MD we want + # to skip real constraints but include special "constraints" + # Like Hookean. + for constraint in self.constraints: + if md and hasattr(constraint, 'redistribute_forces_md'): + constraint.redistribute_forces_md(self, forces) + if not md or hasattr(constraint, 'adjust_potential_energy'): + constraint.adjust_forces(self, forces) + return forces + + # Informs calculators (e.g. Asap) that ideal gas contribution is added here. + _ase_handles_dynamic_stress = True + + def get_stress(self, voigt=True, apply_constraint=True, include_ideal_gas=False): + """Calculate stress tensor. + + Returns an array of the six independent components of the + symmetric stress tensor, in the traditional Voigt order + (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt + order. + + The ideal gas contribution to the stresses is added if the + atoms have momenta and ``include_ideal_gas`` is set to True. + """ + + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + + stress = self._calc.get_stress(self) + shape = stress.shape + + if shape == (3, 3): + # Convert to the Voigt form before possibly applying + # constraints and adding the dynamic part of the stress + # (the "ideal gas contribution"). + stress = full_3x3_to_voigt_6_stress(stress) + else: + assert shape == (6,) + + if apply_constraint: + for constraint in self.constraints: + if hasattr(constraint, 'adjust_stress'): + constraint.adjust_stress(self, stress) + + # Add ideal gas contribution, if applicable + if include_ideal_gas and self.has('momenta'): + stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]]) + p = self.get_momenta() + masses = self.get_masses() + invmass = 1.0 / masses + invvol = 1.0 / self.get_volume() + for alpha in range(3): + for beta in range(alpha, 3): + stress[stresscomp[alpha, beta]] -= ( + p[:, alpha] * p[:, beta] * invmass + ).sum() * invvol + + if voigt: + return stress + else: + return voigt_6_to_full_3x3_stress(stress) + + def get_stresses(self, include_ideal_gas=False, voigt=True): + """Calculate the stress-tensor of all the atoms. + + Only available with calculators supporting per-atom energies and + stresses (e.g. classical potentials). Even for such calculators + there is a certain arbitrariness in defining per-atom stresses. + + The ideal gas contribution to the stresses is added if the + atoms have momenta and ``include_ideal_gas`` is set to True. + """ + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + stresses = self._calc.get_stresses(self) + + # make sure `stresses` are in voigt form + if np.shape(stresses)[1:] == (3, 3): + stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses] + stresses = np.array(stresses_voigt) + + # REMARK: The ideal gas contribution is intensive, i.e., the volume + # is divided out. We currently don't check if `stresses` are intensive + # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`. + # It might be good to check this here, but adds computational overhead. + + if include_ideal_gas and self.has('momenta'): + stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]]) + if hasattr(self._calc, 'get_atomic_volumes'): + invvol = 1.0 / self._calc.get_atomic_volumes() + else: + invvol = self.get_global_number_of_atoms() / self.get_volume() + p = self.get_momenta() + invmass = 1.0 / self.get_masses() + for alpha in range(3): + for beta in range(alpha, 3): + stresses[:, stresscomp[alpha, beta]] -= ( + p[:, alpha] * p[:, beta] * invmass * invvol + ) + if voigt: + return stresses + else: + stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses] + return np.array(stresses_3x3) + + def get_dipole_moment(self): + """Calculate the electric dipole moment for the atoms object. + + Only available for calculators which has a get_dipole_moment() + method.""" + + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + return self._calc.get_dipole_moment(self) + + def copy(self): + """Return a copy.""" + atoms = self.__class__( + cell=self.cell, pbc=self.pbc, info=self.info, celldisp=self._celldisp.copy() + ) + + atoms.arrays = {} + for name, a in self.arrays.items(): + atoms.arrays[name] = a.copy() + atoms.constraints = copy.deepcopy(self.constraints) + return atoms + + def todict(self): + """For basic JSON (non-database) support.""" + d = dict(self.arrays) + d['cell'] = np.asarray(self.cell) + d['pbc'] = self.pbc + if self._celldisp.any(): + d['celldisp'] = self._celldisp + if self.constraints: + d['constraints'] = self.constraints + if self.info: + d['info'] = self.info + # Calculator... trouble. + return d + + @classmethod + def fromdict(cls, dct): + """Rebuild atoms object from dictionary representation (todict).""" + dct = dct.copy() + kw = {} + for name in ['numbers', 'positions', 'cell', 'pbc']: + kw[name] = dct.pop(name) + + constraints = dct.pop('constraints', None) + if constraints: + from ase.constraints import dict2constraint + + constraints = [dict2constraint(d) for d in constraints] + + info = dct.pop('info', None) + + atoms = cls( + constraint=constraints, celldisp=dct.pop('celldisp', None), info=info, **kw + ) + natoms = len(atoms) + + # Some arrays are named differently from the atoms __init__ keywords. + # Also, there may be custom arrays. Hence we set them directly: + for name, arr in dct.items(): + assert len(arr) == natoms, name + assert isinstance(arr, np.ndarray) + atoms.arrays[name] = arr + return atoms + + def __len__(self): + return len(self.arrays['positions']) + + def get_number_of_atoms(self): + """Deprecated, please do not use. + + You probably want len(atoms). Or if your atoms are distributed, + use (and see) get_global_number_of_atoms().""" + import warnings + + warnings.warn( + 'Use get_global_number_of_atoms() instead', np.VisibleDeprecationWarning + ) + return len(self) + + def get_global_number_of_atoms(self): + """Returns the global number of atoms in a distributed-atoms parallel + simulation. + + DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING! + + Equivalent to len(atoms) in the standard ASE Atoms class. You should + normally use len(atoms) instead. This function's only purpose is to + make compatibility between ASE and Asap easier to maintain by having a + few places in ASE use this function instead. It is typically only + when counting the global number of degrees of freedom or in similar + situations. + """ + return len(self) + + def __repr__(self): + tokens = [] + + N = len(self) + if N <= 60: + symbols = self.get_chemical_formula('reduce') + else: + symbols = self.get_chemical_formula('hill') + tokens.append(f"symbols='{symbols}'") + + if self.pbc.any() and not self.pbc.all(): + tokens.append(f'pbc={self.pbc.tolist()}') + else: + tokens.append(f'pbc={self.pbc[0]}') + + cell = self.cell + if cell: + if cell.orthorhombic: + cell = cell.lengths().tolist() + else: + cell = cell.tolist() + tokens.append(f'cell={cell}') + + for name in sorted(self.arrays): + if name in ['numbers', 'positions']: + continue + tokens.append(f'{name}=...') + + if self.constraints: + if len(self.constraints) == 1: + constraint = self.constraints[0] + else: + constraint = self.constraints + tokens.append(f'constraint={repr(constraint)}') + + if self._calc is not None: + tokens.append(f'calculator={self._calc.__class__.__name__}(...)') + + return '{0}({1})'.format(self.__class__.__name__, ', '.join(tokens)) + + def __add__(self, other): + atoms = self.copy() + atoms += other + return atoms + + def extend(self, other): + """Extend atoms object by appending atoms from *other*.""" + if isinstance(other, Atom): + other = self.__class__([other]) + + n1 = len(self) + n2 = len(other) + + for name, a1 in self.arrays.items(): + a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype) + a[:n1] = a1 + if name == 'masses': + a2 = other.get_masses() + else: + a2 = other.arrays.get(name) + if a2 is not None: + a[n1:] = a2 + self.arrays[name] = a + + for name, a2 in other.arrays.items(): + if name in self.arrays: + continue + a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype) + a[n1:] = a2 + if name == 'masses': + a[:n1] = self.get_masses()[:n1] + else: + a[:n1] = 0 + + self.set_array(name, a) + + def __iadd__(self, other): + self.extend(other) + return self + + def append(self, atom): + """Append atom to end.""" + self.extend(self.__class__([atom])) + + def __iter__(self): + for i in range(len(self)): + yield self[i] + + def __getitem__(self, i): + """Return a subset of the atoms. + + i -- scalar integer, list of integers, or slice object + describing which atoms to return. + + If i is a scalar, return an Atom object. If i is a list or a + slice, return an Atoms object with the same cell, pbc, and + other associated info as the original Atoms object. The + indices of the constraints will be shuffled so that they match + the indexing in the subset returned. + + """ + + if isinstance(i, numbers.Integral): + natoms = len(self) + if i < -natoms or i >= natoms: + raise IndexError('Index out of range.') + + return Atom(atoms=self, index=i) + elif not isinstance(i, slice): + i = np.array(i) + # if i is a mask + if i.dtype == bool: + if len(i) != len(self): + raise IndexError( + f'Length of mask {len(i)} must equal ' + f'number of atoms {len(self)}' + ) + i = np.arange(len(self))[i] + + import copy + + conadd = [] + # Constraints need to be deepcopied, but only the relevant ones. + for con in copy.deepcopy(self.constraints): + try: + con.index_shuffle(self, i) + except (IndexError, NotImplementedError): + pass + else: + conadd.append(con) + + atoms = self.__class__( + cell=self.cell, + pbc=self.pbc, + info=self.info, + # should be communicated to the slice as well + celldisp=self._celldisp, + ) + # TODO: Do we need to shuffle indices in adsorbate_info too? + + atoms.arrays = {} + for name, a in self.arrays.items(): + atoms.arrays[name] = a[i].copy() + + atoms.constraints = conadd + return atoms + + def __delitem__(self, i): + from ase.constraints import FixAtoms + + for c in self._constraints: + if not isinstance(c, FixAtoms): + raise RuntimeError( + 'Remove constraint using set_constraint() ' 'before deleting atoms.' + ) + + if isinstance(i, list) and len(i) > 0: + # Make sure a list of booleans will work correctly and not be + # interpreted at 0 and 1 indices. + i = np.array(i) + + if len(self._constraints) > 0: + n = len(self) + i = np.arange(n)[i] + if isinstance(i, int): + i = [i] + constraints = [] + for c in self._constraints: + c = c.delete_atoms(i, n) + if c is not None: + constraints.append(c) + self.constraints = constraints + + mask = np.ones(len(self), bool) + mask[i] = False + for name, a in self.arrays.items(): + self.arrays[name] = a[mask] + + def pop(self, i=-1): + """Remove and return atom at index *i* (default last).""" + atom = self[i] + atom.cut_reference_to_atoms() + del self[i] + return atom + + def __imul__(self, m): + """In-place repeat of atoms.""" + if isinstance(m, int): + m = (m, m, m) + + for x, vec in zip(m, self.cell): + if x != 1 and not vec.any(): + raise ValueError('Cannot repeat along undefined lattice ' 'vector') + + M = np.product(m) + n = len(self) + + for name, a in self.arrays.items(): + self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1)) + + positions = self.arrays['positions'] + i0 = 0 + for m0 in range(m[0]): + for m1 in range(m[1]): + for m2 in range(m[2]): + i1 = i0 + n + positions[i0:i1] += np.dot((m0, m1, m2), self.cell) + i0 = i1 + + if self.constraints is not None: + self.constraints = [c.repeat(m, n) for c in self.constraints] + + self.cell = np.array([m[c] * self.cell[c] for c in range(3)]) + + return self + + def repeat(self, rep): + """Create new repeated atoms object. + + The *rep* argument should be a sequence of three positive + integers like *(2,3,1)* or a single integer (*r*) equivalent + to *(r,r,r)*.""" + + atoms = self.copy() + atoms *= rep + return atoms + + def __mul__(self, rep): + return self.repeat(rep) + + def translate(self, displacement): + """Translate atomic positions. + + The displacement argument can be a float an xyz vector or an + nx3 array (where n is the number of atoms).""" + + self.arrays['positions'] += np.array(displacement) + + def center(self, vacuum=None, axis=(0, 1, 2), about=None): + """Center atoms in unit cell. + + Centers the atoms in the unit cell, so there is the same + amount of vacuum on all sides. + + vacuum: float (default: None) + If specified adjust the amount of vacuum when centering. + If vacuum=10.0 there will thus be 10 Angstrom of vacuum + on each side. + axis: int or sequence of ints + Axis or axes to act on. Default: Act on all axes. + about: float or array (default: None) + If specified, center the atoms about . + I.e., about=(0., 0., 0.) (or just "about=0.", interpreted + identically), to center about the origin. + """ + + # Find the orientations of the faces of the unit cell + cell = self.cell.complete() + dirs = np.zeros_like(cell) + + lengths = cell.lengths() + for i in range(3): + dirs[i] = np.cross(cell[i - 1], cell[i - 2]) + dirs[i] /= np.linalg.norm(dirs[i]) + if dirs[i] @ cell[i] < 0.0: + dirs[i] *= -1 + + if isinstance(axis, int): + axes = (axis,) + else: + axes = axis + + # Now, decide how much each basis vector should be made longer + pos = self.positions + longer = np.zeros(3) + shift = np.zeros(3) + for i in axes: + if len(pos): + scalarprod = pos @ dirs[i] + p0 = scalarprod.min() + p1 = scalarprod.max() + else: + p0 = 0 + p1 = 0 + height = cell[i] @ dirs[i] + if vacuum is not None: + lng = (p1 - p0 + 2 * vacuum) - height + else: + lng = 0.0 # Do not change unit cell size! + top = lng + height - p1 + shf = 0.5 * (top - p0) + cosphi = cell[i] @ dirs[i] / lengths[i] + longer[i] = lng / cosphi + shift[i] = shf / cosphi + + # Now, do it! + translation = np.zeros(3) + for i in axes: + nowlen = lengths[i] + if vacuum is not None: + self.cell[i] = cell[i] * (1 + longer[i] / nowlen) + translation += shift[i] * cell[i] / nowlen + + # We calculated translations using the completed cell, + # so directions without cell vectors will have been centered + # along a "fake" vector of length 1. + # Therefore, we adjust by -0.5: + if not any(self.cell[i]): + translation[i] -= 0.5 + + # Optionally, translate to center about a point in space. + if about is not None: + for vector in self.cell: + translation -= vector / 2.0 + translation += about + + self.positions += translation + + def get_center_of_mass(self, scaled=False): + """Get the center of mass. + + If scaled=True the center of mass in scaled coordinates + is returned.""" + masses = self.get_masses() + com = masses @ self.positions / masses.sum() + if scaled: + return self.cell.scaled_positions(com) + else: + return com + + def set_center_of_mass(self, com, scaled=False): + """Set the center of mass. + + If scaled=True the center of mass is expected in scaled coordinates. + Constraints are considered for scaled=False. + """ + old_com = self.get_center_of_mass(scaled=scaled) + difference = old_com - com + if scaled: + self.set_scaled_positions(self.get_scaled_positions() + difference) + else: + self.set_positions(self.get_positions() + difference) + + def get_moments_of_inertia(self, vectors=False): + """Get the moments of inertia along the principal axes. + + The three principal moments of inertia are computed from the + eigenvalues of the symmetric inertial tensor. Periodic boundary + conditions are ignored. Units of the moments of inertia are + amu*angstrom**2. + """ + com = self.get_center_of_mass() + positions = self.get_positions() + positions -= com # translate center of mass to origin + masses = self.get_masses() + + # Initialize elements of the inertial tensor + I11 = I22 = I33 = I12 = I13 = I23 = 0.0 + for i in range(len(self)): + x, y, z = positions[i] + m = masses[i] + + I11 += m * (y**2 + z**2) + I22 += m * (x**2 + z**2) + I33 += m * (x**2 + y**2) + I12 += -m * x * y + I13 += -m * x * z + I23 += -m * y * z + + I = np.array([[I11, I12, I13], [I12, I22, I23], [I13, I23, I33]]) + + evals, evecs = np.linalg.eigh(I) + if vectors: + return evals, evecs.transpose() + else: + return evals + + def get_angular_momentum(self): + """Get total angular momentum with respect to the center of mass.""" + com = self.get_center_of_mass() + positions = self.get_positions() + positions -= com # translate center of mass to origin + return np.cross(positions, self.get_momenta()).sum(0) + + def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False): + """Rotate atoms based on a vector and an angle, or two vectors. + + Parameters: + + a = None: + Angle that the atoms is rotated around the vector 'v'. 'a' + can also be a vector and then 'a' is rotated + into 'v'. + + v: + Vector to rotate the atoms around. Vectors can be given as + strings: 'x', '-x', 'y', ... . + + center = (0, 0, 0): + The center is kept fixed under the rotation. Use 'COM' to fix + the center of mass, 'COP' to fix the center of positions or + 'COU' to fix the center of cell. + + rotate_cell = False: + If true the cell is also rotated. + + Examples: + + Rotate 90 degrees around the z-axis, so that the x-axis is + rotated into the y-axis: + + >>> atoms = Atoms() + >>> atoms.rotate(90, 'z') + >>> atoms.rotate(90, (0, 0, 1)) + >>> atoms.rotate(-90, '-z') + >>> atoms.rotate('x', 'y') + >>> atoms.rotate((1, 0, 0), (0, 1, 0)) + """ + + if not isinstance(a, numbers.Real): + a, v = v, a + + norm = np.linalg.norm + v = string2vector(v) + + normv = norm(v) + + if normv == 0.0: + raise ZeroDivisionError('Cannot rotate: norm(v) == 0') + + if isinstance(a, numbers.Real): + a *= pi / 180 + v /= normv + c = cos(a) + s = sin(a) + else: + v2 = string2vector(a) + v /= normv + normv2 = np.linalg.norm(v2) + if normv2 == 0: + raise ZeroDivisionError('Cannot rotate: norm(a) == 0') + v2 /= norm(v2) + c = np.dot(v, v2) + v = np.cross(v, v2) + s = norm(v) + # In case *v* and *a* are parallel, np.cross(v, v2) vanish + # and can't be used as a rotation axis. However, in this + # case any rotation axis perpendicular to v2 will do. + eps = 1e-7 + if s < eps: + v = np.cross((0, 0, 1), v2) + if norm(v) < eps: + v = np.cross((1, 0, 0), v2) + assert norm(v) >= eps + elif s > 0: + v /= s + + center = self._centering_as_array(center) + + p = self.arrays['positions'] - center + self.arrays['positions'][:] = ( + c * p - np.cross(p, s * v) + np.outer(np.dot(p, v), (1.0 - c) * v) + center + ) + if rotate_cell: + rotcell = self.get_cell() + rotcell[:] = ( + c * rotcell + - np.cross(rotcell, s * v) + + np.outer(np.dot(rotcell, v), (1.0 - c) * v) + ) + self.set_cell(rotcell) + + def _centering_as_array(self, center): + if isinstance(center, str): + if center.lower() == 'com': + center = self.get_center_of_mass() + elif center.lower() == 'cop': + center = self.get_positions().mean(axis=0) + elif center.lower() == 'cou': + center = self.get_cell().sum(axis=0) / 2 + else: + raise ValueError('Cannot interpret center') + else: + center = np.array(center, float) + return center + + def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)): + """Rotate atoms via Euler angles (in degrees). + + See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation. + + Parameters: + + center : + The point to rotate about. A sequence of length 3 with the + coordinates, or 'COM' to select the center of mass, 'COP' to + select center of positions or 'COU' to select center of cell. + phi : + The 1st rotation angle around the z axis. + theta : + Rotation around the x axis. + psi : + 2nd rotation around the z axis. + + """ + center = self._centering_as_array(center) + + phi *= pi / 180 + theta *= pi / 180 + psi *= pi / 180 + + # First move the molecule to the origin In contrast to MATLAB, + # numpy broadcasts the smaller array to the larger row-wise, + # so there is no need to play with the Kronecker product. + rcoords = self.positions - center + # First Euler rotation about z in matrix form + D = np.array( + ((cos(phi), sin(phi), 0.0), (-sin(phi), cos(phi), 0.0), (0.0, 0.0, 1.0)) + ) + # Second Euler rotation about x: + C = np.array( + ( + (1.0, 0.0, 0.0), + (0.0, cos(theta), sin(theta)), + (0.0, -sin(theta), cos(theta)), + ) + ) + # Third Euler rotation, 2nd rotation about z: + B = np.array( + ((cos(psi), sin(psi), 0.0), (-sin(psi), cos(psi), 0.0), (0.0, 0.0, 1.0)) + ) + # Total Euler rotation + A = np.dot(B, np.dot(C, D)) + # Do the rotation + rcoords = np.dot(A, np.transpose(rcoords)) + # Move back to the rotation point + self.positions = np.transpose(rcoords) + center + + def get_dihedral(self, a0, a1, a2, a3, mic=False): + """Calculate dihedral angle. + + Calculate dihedral angle (in degrees) between the vectors a0->a1 + and a2->a3. + + Use mic=True to use the Minimum Image Convention and calculate the + angle across periodic boundaries. + """ + return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0] + + def get_dihedrals(self, indices, mic=False): + """Calculate dihedral angles. + + Calculate dihedral angles (in degrees) between the list of vectors + a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices. + + Use mic=True to use the Minimum Image Convention and calculate the + angles across periodic boundaries. + """ + indices = np.array(indices) + assert indices.shape[1] == 4 + + a0s = self.positions[indices[:, 0]] + a1s = self.positions[indices[:, 1]] + a2s = self.positions[indices[:, 2]] + a3s = self.positions[indices[:, 3]] + + # vectors 0->1, 1->2, 2->3 + v0 = a1s - a0s + v1 = a2s - a1s + v2 = a3s - a2s + + cell = None + pbc = None + + if mic: + cell = self.cell + pbc = self.pbc + + return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc) + + def _masked_rotate(self, center, axis, diff, mask): + # do rotation of subgroup by copying it to temporary atoms object + # and then rotating that + # + # recursive object definition might not be the most elegant thing, + # more generally useful might be a rotation function with a mask? + group = self.__class__() + for i in range(len(self)): + if mask[i]: + group += self[i] + group.translate(-center) + group.rotate(diff * 180 / pi, axis) + group.translate(center) + # set positions in original atoms object + j = 0 + for i in range(len(self)): + if mask[i]: + self.positions[i] = group[j].position + j += 1 + + def set_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None): + """Set the dihedral angle (degrees) between vectors a1->a2 and + a3->a4 by changing the atom indexed by a4. + + If mask is not None, all the atoms described in mask + (read: the entire subgroup) are moved. Alternatively to the mask, + the indices of the atoms to be rotated can be supplied. If both + *mask* and *indices* are given, *indices* overwrites *mask*. + + **Important**: If *mask* or *indices* is given and does not contain + *a4*, *a4* will NOT be moved. In most cases you therefore want + to include *a4* in *mask*/*indices*. + + Example: the following defines a very crude + ethane-like molecule and twists one half of it by 30 degrees. + + >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0], + ... [1, 0, 0], [2, 1, 0], [2, -1, 0]]) + >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1]) + """ + + angle *= pi / 180 + + # if not provided, set mask to the last atom in the + # dihedral description + if mask is None and indices is None: + mask = np.zeros(len(self)) + mask[a4] = 1 + elif indices is not None: + mask = [index in indices for index in range(len(self))] + + # compute necessary in dihedral change, from current value + current = self.get_dihedral(a1, a2, a3, a4) * pi / 180 + diff = angle - current + axis = self.positions[a3] - self.positions[a2] + center = self.positions[a3] + self._masked_rotate(center, axis, diff, mask) + + def rotate_dihedral(self, a1, a2, a3, a4, angle=None, mask=None, indices=None): + """Rotate dihedral angle. + + Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a + predefined dihedral angle, starting from its current configuration. + """ + start = self.get_dihedral(a1, a2, a3, a4) + self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices) + + def get_angle(self, a1, a2, a3, mic=False): + """Get angle formed by three atoms. + + Calculate angle in degrees between the vectors a2->a1 and + a2->a3. + + Use mic=True to use the Minimum Image Convention and calculate the + angle across periodic boundaries. + """ + return self.get_angles([[a1, a2, a3]], mic=mic)[0] + + def get_angles(self, indices, mic=False): + """Get angle formed by three atoms for multiple groupings. + + Calculate angle in degrees between vectors between atoms a2->a1 + and a2->a3, where a1, a2, and a3 are in each row of indices. + + Use mic=True to use the Minimum Image Convention and calculate + the angle across periodic boundaries. + """ + indices = np.array(indices) + assert indices.shape[1] == 3 + + a1s = self.positions[indices[:, 0]] + a2s = self.positions[indices[:, 1]] + a3s = self.positions[indices[:, 2]] + + v12 = a1s - a2s + v32 = a3s - a2s + + cell = None + pbc = None + + if mic: + cell = self.cell + pbc = self.pbc + + return get_angles(v12, v32, cell=cell, pbc=pbc) + + def set_angle( + self, a1, a2=None, a3=None, angle=None, mask=None, indices=None, add=False + ): + """Set angle (in degrees) formed by three atoms. + + Sets the angle between vectors *a2*->*a1* and *a2*->*a3*. + + If *add* is `True`, the angle will be changed by the value given. + + Same usage as in :meth:`ase.Atoms.set_dihedral`. + If *mask* and *indices* + are given, *indices* overwrites *mask*. If *mask* and *indices* + are not set, only *a3* is moved.""" + + if any(a is None for a in [a2, a3, angle]): + raise ValueError('a2, a3, and angle must not be None') + + # If not provided, set mask to the last atom in the angle description + if mask is None and indices is None: + mask = np.zeros(len(self)) + mask[a3] = 1 + elif indices is not None: + mask = [index in indices for index in range(len(self))] + + if add: + diff = angle + else: + # Compute necessary in angle change, from current value + diff = angle - self.get_angle(a1, a2, a3) + + diff *= pi / 180 + # Do rotation of subgroup by copying it to temporary atoms object and + # then rotating that + v10 = self.positions[a1] - self.positions[a2] + v12 = self.positions[a3] - self.positions[a2] + v10 /= np.linalg.norm(v10) + v12 /= np.linalg.norm(v12) + axis = np.cross(v10, v12) + center = self.positions[a2] + self._masked_rotate(center, axis, diff, mask) + + def rattle(self, stdev=0.001, seed=None, rng=None): + """Randomly displace atoms. + + This method adds random displacements to the atomic positions, + taking a possible constraint into account. The random numbers are + drawn from a normal distribution of standard deviation stdev. + + For a parallel calculation, it is important to use the same + seed on all processors!""" + + if seed is not None and rng is not None: + raise ValueError('Please do not provide both seed and rng.') + + if rng is None: + if seed is None: + seed = 42 + rng = np.random.RandomState(seed) + positions = self.arrays['positions'] + self.set_positions(positions + rng.normal(scale=stdev, size=positions.shape)) + + def get_distance(self, a0, a1, mic=False, vector=False): + """Return distance between two atoms. + + Use mic=True to use the Minimum Image Convention. + vector=True gives the distance vector (from a0 to a1). + """ + return self.get_distances(a0, [a1], mic=mic, vector=vector)[0] + + def get_distances(self, a, indices, mic=False, vector=False): + """Return distances of atom No.i with a list of atoms. + + Use mic=True to use the Minimum Image Convention. + vector=True gives the distance vector (from a to self[indices]). + """ + R = self.arrays['positions'] + p1 = [R[a]] + p2 = R[indices] + + cell = None + pbc = None + + if mic: + cell = self.cell + pbc = self.pbc + + D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc) + + if vector: + D.shape = (-1, 3) + return D + else: + D_len.shape = (-1,) + return D_len + + def get_all_distances(self, mic=False, vector=False): + """Return distances of all of the atoms with all of the atoms. + + Use mic=True to use the Minimum Image Convention. + """ + R = self.arrays['positions'] + + cell = None + pbc = None + + if mic: + cell = self.cell + pbc = self.pbc + + D, D_len = get_distances(R, cell=cell, pbc=pbc) + + if vector: + return D + else: + return D_len + + def set_distance( + self, + a0, + a1, + distance, + fix=0.5, + mic=False, + mask=None, + indices=None, + add=False, + factor=False, + ): + """Set the distance between two atoms. + + Set the distance between atoms *a0* and *a1* to *distance*. + By default, the center of the two atoms will be fixed. Use + *fix=0* to fix the first atom, *fix=1* to fix the second + atom and *fix=0.5* (default) to fix the center of the bond. + + If *mask* or *indices* are set (*mask* overwrites *indices*), + only the atoms defined there are moved + (see :meth:`ase.Atoms.set_dihedral`). + + When *add* is true, the distance is changed by the value given. + In combination + with *factor* True, the value given is a factor scaling the distance. + + It is assumed that the atoms in *mask*/*indices* move together + with *a1*. If *fix=1*, only *a0* will therefore be moved.""" + + if a0 % len(self) == a1 % len(self): + raise ValueError('a0 and a1 must not be the same') + + if add: + oldDist = self.get_distance(a0, a1, mic=mic) + if factor: + newDist = oldDist * distance + else: + newDist = oldDist + distance + self.set_distance( + a0, + a1, + newDist, + fix=fix, + mic=mic, + mask=mask, + indices=indices, + add=False, + factor=False, + ) + return + + R = self.arrays['positions'] + D = np.array([R[a1] - R[a0]]) + + if mic: + D, D_len = find_mic(D, self.cell, self.pbc) + else: + D_len = np.array([np.sqrt((D**2).sum())]) + x = 1.0 - distance / D_len[0] + + if mask is None and indices is None: + indices = [a0, a1] + elif mask: + indices = [i for i in range(len(self)) if mask[i]] + + for i in indices: + if i == a0: + R[a0] += (x * fix) * D[0] + else: + R[i] -= (x * (1.0 - fix)) * D[0] + + def get_scaled_positions(self, wrap=True): + """Get positions relative to unit cell. + + If wrap is True, atoms outside the unit cell will be wrapped into + the cell in those directions with periodic boundary conditions + so that the scaled coordinates are between zero and one. + + If any cell vectors are zero, the corresponding coordinates + are evaluated as if the cell were completed using + ``cell.complete()``. This means coordinates will be Cartesian + as long as the non-zero cell vectors span a Cartesian axis or + plane.""" + + fractional = self.cell.scaled_positions(self.positions) + + if wrap: + for i, periodic in enumerate(self.pbc): + if periodic: + # Yes, we need to do it twice. + # See the scaled_positions.py test. + fractional[:, i] %= 1.0 + fractional[:, i] %= 1.0 + + return fractional + + def set_scaled_positions(self, scaled): + """Set positions relative to unit cell.""" + self.positions[:] = self.cell.cartesian_positions(scaled) + + def wrap(self, **wrap_kw): + """Wrap positions to unit cell. + + Parameters: + + wrap_kw: (keyword=value) pairs + optional keywords `pbc`, `center`, `pretty_translation`, `eps`, + see :func:`ase.geometry.wrap_positions` + """ + + if 'pbc' not in wrap_kw: + wrap_kw['pbc'] = self.pbc + + self.positions[:] = self.get_positions(wrap=True, **wrap_kw) + + def get_temperature(self): + """Get the temperature in Kelvin.""" + dof = len(self) * 3 + for constraint in self._constraints: + dof -= constraint.get_removed_dof(self) + ekin = self.get_kinetic_energy() + return 2 * ekin / (dof * units.kB) + + def __eq__(self, other): + """Check for identity of two atoms objects. + + Identity means: same positions, atomic numbers, unit cell and + periodic boundary conditions.""" + if not isinstance(other, Atoms): + return False + a = self.arrays + b = other.arrays + return ( + len(self) == len(other) + and (a['positions'] == b['positions']).all() + and (a['numbers'] == b['numbers']).all() + and (self.cell == other.cell).all() + and (self.pbc == other.pbc).all() + ) + + def __ne__(self, other): + """Check if two atoms objects are not equal. + + Any differences in positions, atomic numbers, unit cell or + periodic boundary condtions make atoms objects not equal. + """ + eq = self.__eq__(other) + if eq is NotImplemented: + return eq + else: + return not eq + + # @deprecated('Please use atoms.cell.volume') + # We kind of want to deprecate this, but the ValueError behaviour + # might be desirable. Should we do this? + def get_volume(self): + """Get volume of unit cell.""" + if self.cell.rank != 3: + raise ValueError( + f'You have {self.cell.rank} lattice vectors: volume not defined' + ) + return self.cell.volume + + def _get_positions(self): + """Return reference to positions-array for in-place manipulations.""" + return self.arrays['positions'] + + def _set_positions(self, pos): + """Set positions directly, bypassing constraints.""" + self.arrays['positions'][:] = pos + + positions = property( + _get_positions, + _set_positions, + doc='Attribute for direct ' + 'manipulation of the positions.', + ) + + def _get_atomic_numbers(self): + """Return reference to atomic numbers for in-place + manipulations.""" + return self.arrays['numbers'] + + numbers = property( + _get_atomic_numbers, + set_atomic_numbers, + doc='Attribute for direct ' + 'manipulation of the atomic numbers.', + ) + + @property + def cell(self): + """The :class:`ase.cell.Cell` for direct manipulation.""" + return self._cellobj + + @cell.setter + def cell(self, cell): + cell = Cell.ascell(cell) + self._cellobj[:] = cell + + def write(self, filename, format=None, **kwargs): + """Write atoms object to a file. + + see ase.io.write for formats. + kwargs are passed to ase.io.write. + """ + from ase.io import write + + write(filename, self, format, **kwargs) + + def iterimages(self): + yield self + + def edit(self): + """Modify atoms interactively through ASE's GUI viewer. + + Conflicts leading to undesirable behaviour might arise + when matplotlib has been pre-imported with certain + incompatible backends and while trying to use the + plot feature inside the interactive GUI. To circumvent, + please set matplotlib.use('gtk') before calling this + method. + """ + from ase.gui.gui import GUI + from ase.gui.images import Images + + images = Images([self]) + gui = GUI(images) + gui.run() + + +def string2vector(v): + if isinstance(v, str): + if v[0] == '-': + return -string2vector(v[1:]) + w = np.zeros(3) + w['xyz'.index(v)] = 1.0 + return w + return np.array(v, float) + + +def default(data, dflt): + """Helper function for setting default values.""" + if data is None: + return None + elif isinstance(data, (list, tuple)): + newdata = [] + allnone = True + for x in data: + if x is None: + newdata.append(dflt) + else: + newdata.append(x) + allnone = False + if allnone: + return None + return newdata + else: + return data + + +# ? How generic (usable for any CG model) vs. Martini-specific do we want to be? +class ParticlesState(Entity): + """ + A base section to define individual coarse-grained (CG) particle information. + """ + + # ? What do we want to qualify as type identifier? What safety checks do we need? + particle_type = Quantity( + type=str, + description=""" + Symbol(s) describing the CG particle type. Currently, entrie particle label is + used for type definition. + """, + ) + + # ? Do we want to reflect the Martini size nomenclature and include bead volume/bead mass? + # particle_size = Quantity( + # type=np.float64, + # description=""" + # Particle size, determining the number of non-hydrogen atoms represented by the + # particle. Currently, possible values are 0.47 nm (regular, default), + # 0.43/0.41 nm (small), and 0.34 nm (tiny). + # """, + # ) + + # particle_mass = Quantity( + # type=np.float64, + # description=""" + # Particle size, determining the number of non-hydrogen atoms represented by the + # particle. Currently, possible values are 72 amu (regular, default), 54/45 amu + # (small), and 36 amu (tiny). + # """, + # ) + + charge = Quantity( + type=np.int32, + default=0, + description=""" + Charge of the particle. Neutral = 0. Can be any positive integer (+1, +2...) + for cations or any negative integer (-1, -2...) for anions. + """, + a_eln=ELNAnnotation(component='NumberEditQuantity'), + ) + + def resolve_particle_type(self, logger: 'BoundLogger') -> Optional[str]: + """ + Checks if any value is passed as particle label. Converts to string to be used as + type identifier for the CG particle. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[str]): The resolved `particle type`. + """ + if self.particle_type is not None and self.particle_type.isascii(): + try: + return str(self.particle_type) + except TypeError: + logger.error('The parsed `particle type` can not be read.') + return None + + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) + + # Get particle_type as string, if possible. + if not isinstance(self.particle_type, str): + self.particle_type = self.resolve_particle_type(logger=logger) From beb7a8b7d9412cc18d11c8d2dea8a9dec8b8713e Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Fri, 11 Oct 2024 15:28:21 +0200 Subject: [PATCH 05/16] working on CG schema --- .../schema_packages/model_system.py | 162 +++ .../schema_packages/particles_state.py | 1077 ++++++----------- 2 files changed, 499 insertions(+), 740 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 0555c432..4671bf56 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -27,6 +27,7 @@ from structlog.stdlib import BoundLogger from nomad_simulations.schema_packages.atoms_state import AtomsState +from nomad_simulations.schema_packages.particles_state import ParticlesState from nomad_simulations.schema_packages.utils import ( get_sibling_section, is_not_representative, @@ -492,6 +493,167 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: self.name = self.m_def.name if self.name is None else self.name +class ParticleCell(Cell): + """ + A base section used to specify the atomic cell information of a system. + """ + + particles_state = SubSection(sub_section=ParticlesState.m_def, repeats=True) + + n_particles = Quantity( + type=np.int32, + description=""" + Number of atoms in the atomic cell. + """, + ) + + equivalent_particles = Quantity( + type=np.int32, + shape=['n_atoms'], + description=""" + List of equivalent atoms as defined in `atoms`. If no equivalent atoms are found, + then the list is simply the index of each element, e.g.: + - [0, 1, 2, 3] all four atoms are non-equivalent. + - [0, 0, 0, 3] three equivalent atoms and one non-equivalent. + """, + ) + + def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs): + super().__init__(m_def, m_context, **kwargs) + # Set the name of the section + self.name = self.m_def.name + + def is_equal_cell(self, other) -> bool: + """ + Check if the atomic cell is equal to an`other` atomic cell by comparing the `positions` and + the `AtomsState[*].chemical_symbol`. + Args: + other: The other atomic cell to compare with. + Returns: + bool: True if the atomic cells are equal, False otherwise. + """ + if not isinstance(other, ParticleCell): + return False + + # Compare positions using the parent sections's `__eq__` method + if not super().is_equal_cell(other=other): + return False + + # Check that the `chemical_symbol` of the atoms in `cell_1` match with the ones in `cell_2` + check_positions = self._check_positions( + positions_1=self.positions, positions_2=other.positions + ) + try: + for particle in check_positions: + type_1 = self.particles_state[particle[0]].particle_type + type_2 = other.particles_state[particle[1]].particle_type + if type_1 != type_2: + return False + except Exception: + return False + return True + + # def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]: + # """ + # Get the chemical symbols of the atoms in the atomic cell. These are defined on `atoms_state[*].chemical_symbol`. + + # Args: + # logger (BoundLogger): The logger to log messages. + + # Returns: + # list: The list of chemical symbols of the atoms in the atomic cell. + # """ + # if not self.atoms_state: + # return [] + + # chemical_symbols = [] + # for atom_state in self.atoms_state: + # if not atom_state.chemical_symbol: + # logger.warning('Could not find `AtomsState[*].chemical_symbol`.') + # return [] + # chemical_symbols.append(atom_state.chemical_symbol) + # return chemical_symbols + + # def to_ase_atoms(self, logger: 'BoundLogger') -> Optional[ase.Atoms]: + # """ + # Generates an ASE Atoms object with the most basic information from the parsed `AtomicCell` + # section (labels, periodic_boundary_conditions, positions, and lattice_vectors). + + # Args: + # logger (BoundLogger): The logger to log messages. + + # Returns: + # (Optional[ase.Atoms]): The ASE Atoms object with the basic information from the `AtomicCell`. + # """ + # # Initialize ase.Atoms object with labels + # atoms_labels = self.get_chemical_symbols(logger=logger) + # ase_atoms = ase.Atoms(symbols=atoms_labels) + + # # PBC + # if self.periodic_boundary_conditions is None: + # logger.info( + # 'Could not find `AtomicCell.periodic_boundary_conditions`. They will be set to [False, False, False].' + # ) + # self.periodic_boundary_conditions = [False, False, False] + # ase_atoms.set_pbc(pbc=self.periodic_boundary_conditions) + + # # Lattice vectors + # if self.lattice_vectors is not None: + # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + # else: + # logger.info('Could not find `AtomicCell.lattice_vectors`.') + + # # Positions + # if self.positions is not None: + # if len(self.positions) != len(self.atoms_state): + # logger.error( + # 'Length of `AtomicCell.positions` does not coincide with the length of the `AtomicCell.atoms_state`.' + # ) + # return None + # ase_atoms.set_positions( + # newpositions=self.positions.to('angstrom').magnitude + # ) + # else: + # logger.warning('Could not find `AtomicCell.positions`.') + # return None + + # return ase_atoms + + # def from_ase_atoms(self, ase_atoms: ase.Atoms, logger: 'BoundLogger') -> None: + # """ + # Parses the information from an ASE Atoms object to the `AtomicCell` section. + + # Args: + # ase_atoms (ase.Atoms): The ASE Atoms object to parse. + # logger (BoundLogger): The logger to log messages. + # """ + # # `AtomsState[*].chemical_symbol` + # for symbol in ase_atoms.get_chemical_symbols(): + # atom_state = AtomsState(chemical_symbol=symbol) + # self.atoms_state.append(atom_state) + + # # `periodic_boundary_conditions` + # self.periodic_boundary_conditions = ase_atoms.get_pbc() + + # # `lattice_vectors` + # cell = ase_atoms.get_cell() + # self.lattice_vectors = ase.geometry.complete_cell(cell) * ureg('angstrom') + + # # `positions` + # positions = ase_atoms.get_positions() + # if ( + # not positions.tolist() + # ): # ASE assigns a shape=(0, 3) array if no positions are found + # return None + # self.positions = positions * ureg('angstrom') + + def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: + super().normalize(archive, logger) + + # Set the name of the section + self.name = self.m_def.name if self.name is None else self.name + + class Symmetry(ArchiveSection): """ A base section used to specify the symmetry of the `AtomicCell`. diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index 458a00d3..b601e735 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -1,7 +1,10 @@ +import numbers from typing import TYPE_CHECKING, Any, Optional, Union +import ase import numpy as np import pint +from deprecated import deprecated from nomad.datamodel.data import ArchiveSection from nomad.datamodel.metainfo.annotations import ELNAnnotation from nomad.datamodel.metainfo.basesections import Entity @@ -31,31 +34,53 @@ class Particle: Parameters: - types: str (formula) or list of str - Can be a string formula, a list of symbols or a list of - Atom objects. Examples: 'H2O', 'COPt12', ['H', 'H', 'O'], - [Atom('Ne', (x, y, z)), ...]. - positions: list of xyz-positions - Atomic positions. Anything that can be converted to an - ndarray of shape (n, 3) will do: [(x1,y1,z1), (x2,y2,z2), - ...]. - scaled_positions: list of scaled-positions - Like positions, but given in units of the unit cell. - Can not be set at the same time as positions. - numbers: list of int - Atomic numbers (use only one of symbols/numbers). - tags: list of int - Special purpose tags. - momenta: list of xyz-momenta - Momenta for all atoms. + types: str (formula) or list of str. + Default: [‘A’] + typeid: int + (Optional) Id numbers of corresponding particle types. + Default: 0 + type_shapes: str + Store a per-type shape definition for visualization. + A dictionary is stored for each of the NT types, corresponding + to a shape for visualization of that type. + Default: empty masses: list of float - Atomic masses in atomic units. - magmoms: list of float or list of xyz-values - Magnetic moments. Can be either a single value for each atom - for collinear calculations or three numbers for each atom for - non-collinear calculations. + The mass of each particle. + Default: 1.0 charges: list of float Initial atomic charges. + Default: 0.0 + diameter: float + The diameter of each particle. + Default: 1.0 + body: int + The composite body associated with each particle. The value -1 + indicates no body. + Default: -1 + moment_inertia: float + The moment_inertia of each particle (I_xx, I_yy, I_zz). + This inertia tensor is diagonal in the body frame of the particle. + The default value is for point particles. + Default: 0, 0, 0 + positions: float, list of xyz-positions + Particle positions. Needs to be convertible to an + ndarray of shape (N, 3). + Default: 0, 0, 0 + scaled_positions: list of scaled-positions + Like positions, but given in units of the unit cell. + Can not be set at the same time as positions. + Default: 0, 0, 0 + orientation: float + The orientation of each particle. In scalar + vector notation, + this is (r, a_x, a_y, a_z), where the quaternion is q = r + a_xi + a_yj + a_zk. + A unit quaternion has the property: sqrt(r^2 + a_x^2 + a_y^2 + a_z^2) = 1. + Default: 0, 0, 0, 0 + angmom: float + The angular momentum of each particle as a quaternion. + Default: 0, 0, 0, 0 + image: int + The number of times each particle has wrapped around the box (i_x, i_y, i_z). + Default: 0, 0, 0 cell: 3x3 matrix or length 3 or 6 vector Unit cell vectors. Can also be given as just three numbers for orthorhombic cells, or 6 numbers, where @@ -73,25 +98,6 @@ class Particle: Periodic boundary conditions flags. Examples: True, False, 0, 1, (1, 1, 0), (True, False, False). Default value: False. - constraint: constraint object(s) - Used for applying one or more constraints during structure - optimization. - calculator: calculator object - Used to attach a calculator for calculating energies and atomic - forces. - info: dict of key-value pairs - Dictionary of key-value pairs with additional information - about the system. The following keys may be used by ase: - - - spacegroup: Spacegroup instance - - unit_cell: 'conventional' | 'primitive' | int | 3 ints - - adsorbate_info: Information about special adsorption sites - - Items in the info attribute survives copy and slicing and can - be stored in and retrieved from trajectory files given that the - key is a string, the value is JSON-compatible and, if the value is a - user-defined object, its base class is importable. One should - not make any assumptions about the existence of keys. Examples: @@ -118,225 +124,219 @@ class Particle: ... pbc=(1, 0, 0)) """ - ase_objtype = 'atoms' # For JSONability - def __init__( self, - symbols=None, + types=None, positions=None, - numbers=None, - tags=None, - momenta=None, + typeid=None, + type_shapes=None, + moment_inertia=None, masses=None, - magmoms=None, + angmom=None, charges=None, + diameter=None, + body=None, scaled_positions=None, + orientation=None, + image=None, cell=None, pbc=None, celldisp=None, - constraint=None, - calculator=None, - info=None, - velocities=None, ): - self._cellobj = Cell.new() + self._cellobj = self.set_cell() self._pbc = np.zeros(3, bool) - atoms = None - - if hasattr(symbols, 'get_positions'): - atoms = symbols - symbols = None - elif ( - isinstance(symbols, (list, tuple)) - and len(symbols) > 0 - and isinstance(symbols[0], Atom) - ): - # Get data from a list or tuple of Atom objects: - data = [ - [atom.get_raw(name) for atom in symbols] - for name in [ - 'position', - 'number', - 'tag', - 'momentum', - 'mass', - 'magmom', - 'charge', - ] - ] - atoms = self.__class__(None, *data) - symbols = None - - if atoms is not None: - # Get data from another Atoms object: - if scaled_positions is not None: - raise NotImplementedError - if symbols is None and numbers is None: - numbers = atoms.get_atomic_numbers() - if positions is None: - positions = atoms.get_positions() - if tags is None and atoms.has('tags'): - tags = atoms.get_tags() - if momenta is None and atoms.has('momenta'): - momenta = atoms.get_momenta() - if magmoms is None and atoms.has('initial_magmoms'): - magmoms = atoms.get_initial_magnetic_moments() - if masses is None and atoms.has('masses'): - masses = atoms.get_masses() - if charges is None and atoms.has('initial_charges'): - charges = atoms.get_initial_charges() - if cell is None: - cell = atoms.get_cell() - if celldisp is None: - celldisp = atoms.get_celldisp() - if pbc is None: - pbc = atoms.get_pbc() - if constraint is None: - constraint = [c.copy() for c in atoms.constraints] - if calculator is None: - calculator = atoms.calc - if info is None: - info = copy.deepcopy(atoms.info) - - self.arrays = {} - - if symbols is None: - if numbers is None: - if positions is not None: - natoms = len(positions) - elif scaled_positions is not None: - natoms = len(scaled_positions) - else: - natoms = 0 - numbers = np.zeros(natoms, int) - self.new_array('numbers', numbers, int) - else: - if numbers is not None: - raise TypeError('Use only one of "symbols" and "numbers".') - else: - self.new_array('numbers', symbols2numbers(symbols), int) - - if self.numbers.ndim != 1: - raise ValueError('"numbers" must be 1-dimensional.') - - if cell is None: - cell = np.zeros((3, 3)) - self.set_cell(cell) - - if celldisp is None: - celldisp = np.zeros(shape=(3, 1)) - self.set_celldisp(celldisp) - - if positions is None: - if scaled_positions is None: - positions = np.zeros((len(self.arrays['numbers']), 3)) - else: - assert self.cell.rank == 3 - positions = np.dot(scaled_positions, self.cell) - else: - if scaled_positions is not None: - raise TypeError('Use only one of "symbols" and "numbers".') - self.new_array('positions', positions, float, (3,)) - - self.set_constraint(constraint) - self.set_tags(default(tags, 0)) - self.set_masses(default(masses, None)) - self.set_initial_magnetic_moments(default(magmoms, 0.0)) - self.set_initial_charges(default(charges, 0.0)) - if pbc is None: - pbc = False - self.set_pbc(pbc) - self.set_momenta(default(momenta, (0.0, 0.0, 0.0)), apply_constraint=False) - - if velocities is not None: - if momenta is None: - self.set_velocities(velocities) - else: - raise TypeError('Use only one of "momenta" and "velocities".') - - if info is None: - self.info = {} - else: - self.info = dict(info) - - self.calc = calculator - - @property - def symbols(self): - """Get chemical symbols as a :class:`ase.symbols.Symbols` object. - - The object works like ``atoms.numbers`` except its values - are strings. It supports in-place editing.""" - return Symbols(self.numbers) - - @symbols.setter - def symbols(self, obj): - new_symbols = Symbols.fromsymbols(obj) - self.numbers[:] = new_symbols.numbers - - @deprecated(DeprecationWarning('Please use atoms.calc = calc')) - def set_calculator(self, calc=None): - """Attach calculator object. - - Please use the equivalent atoms.calc = calc instead of this - method.""" - self.calc = calc - - @deprecated(DeprecationWarning('Please use atoms.calc')) - def get_calculator(self): - """Get currently attached calculator object. - - Please use the equivalent atoms.calc instead of - atoms.get_calculator().""" - return self.calc - - @property - def calc(self): - """Calculator object.""" - return self._calc - - @calc.setter - def calc(self, calc): - self._calc = calc - if hasattr(calc, 'set_atoms'): - calc.set_atoms(self) - - @calc.deleter # type: ignore - @deprecated(DeprecationWarning('Please use atoms.calc = None')) - def calc(self): - self._calc = None - - @property # type: ignore - @deprecated('Please use atoms.cell.rank instead') - def number_of_lattice_vectors(self): - """Number of (non-zero) lattice vectors.""" - return self.cell.rank - - def set_constraint(self, constraint=None): - """Apply one or more constrains. - - The *constraint* argument must be one constraint object or a - list of constraint objects.""" - if constraint is None: - self._constraints = [] - else: - if isinstance(constraint, list): - self._constraints = constraint - elif isinstance(constraint, tuple): - self._constraints = list(constraint) - else: - self._constraints = [constraint] - - def _get_constraints(self): - return self._constraints - - def _del_constraints(self): - self._constraints = [] - - constraints = property( - _get_constraints, set_constraint, _del_constraints, 'Constraints of the atoms.' - ) + particles = None + + # if hasattr(symbols, 'get_positions'): + # atoms = symbols + # symbols = None + # elif ( + # isinstance(symbols, (list, tuple)) + # and len(symbols) > 0 + # and isinstance(symbols[0], Atom) + # ): + # # Get data from a list or tuple of Atom objects: + # data = [ + # [atom.get_raw(name) for atom in symbols] + # for name in [ + # 'position', + # 'number', + # 'tag', + # 'momentum', + # 'mass', + # 'magmom', + # 'charge', + # ] + # ] + # atoms = self.__class__(None, *data) + # symbols = None + + # if atoms is not None: + # # Get data from another Atoms object: + # if scaled_positions is not None: + # raise NotImplementedError + # if symbols is None and numbers is None: + # numbers = atoms.get_atomic_numbers() + # if positions is None: + # positions = atoms.get_positions() + # if tags is None and atoms.has('tags'): + # tags = atoms.get_tags() + # if momenta is None and atoms.has('momenta'): + # momenta = atoms.get_momenta() + # if magmoms is None and atoms.has('initial_magmoms'): + # magmoms = atoms.get_initial_magnetic_moments() + # if masses is None and atoms.has('masses'): + # masses = atoms.get_masses() + # if charges is None and atoms.has('initial_charges'): + # charges = atoms.get_initial_charges() + # if cell is None: + # cell = atoms.get_cell() + # if celldisp is None: + # celldisp = atoms.get_celldisp() + # if pbc is None: + # pbc = atoms.get_pbc() + + # self.arrays = {} + + # if symbols is None: + # if numbers is None: + # if positions is not None: + # natoms = len(positions) + # elif scaled_positions is not None: + # natoms = len(scaled_positions) + # else: + # natoms = 0 + # numbers = np.zeros(natoms, int) + # self.new_array('numbers', numbers, int) + # else: + # if numbers is not None: + # raise TypeError('Use only one of "symbols" and "numbers".') + # else: + # self.new_array('numbers', symbols2numbers(symbols), int) + + # if self.numbers.ndim != 1: + # raise ValueError('"numbers" must be 1-dimensional.') + + # if cell is None: + # cell = np.zeros((3, 3)) + # self.set_cell(cell) + + # if celldisp is None: + # celldisp = np.zeros(shape=(3, 1)) + # self.set_celldisp(celldisp) + + # if positions is None: + # if scaled_positions is None: + # positions = np.zeros((len(self.arrays['numbers']), 3)) + # else: + # assert self.cell.rank == 3 + # positions = np.dot(scaled_positions, self.cell) + # else: + # if scaled_positions is not None: + # raise TypeError('Use only one of "symbols" and "numbers".') + # self.new_array('positions', positions, float, (3,)) + # self.set_tags(default(tags, 0)) + # self.set_masses(default(masses, None)) + # self.set_initial_magnetic_moments(default(magmoms, 0.0)) + # self.set_initial_charges(default(charges, 0.0)) + # if pbc is None: + # pbc = False + # self.set_pbc(pbc) + # self.set_momenta(default(momenta, (0.0, 0.0, 0.0)), apply_constraint=False) + + # if velocities is not None: + # if momenta is None: + # self.set_velocities(velocities) + # else: + # raise TypeError('Use only one of "momenta" and "velocities".') + + # if info is None: + # self.info = {} + # else: + # self.info = dict(info) + + # self.calc = calculator + + # def set_cell(self, cell): + # if cell is None: + # cell = np.zeros((3, 3)) + + # @property + # def symbols(self): + # """Get chemical symbols as a :class:`ase.symbols.Symbols` object. + + # The object works like ``atoms.numbers`` except its values + # are strings. It supports in-place editing.""" + # return Symbols(self.numbers) + + # @symbols.setter + # def symbols(self, obj): + # new_symbols = Symbols.fromsymbols(obj) + # self.numbers[:] = new_symbols.numbers + + # @deprecated(DeprecationWarning('Please use atoms.calc = calc')) + # def set_calculator(self, calc=None): + # """Attach calculator object. + + # Please use the equivalent atoms.calc = calc instead of this + # method.""" + # self.calc = calc + + # @deprecated(DeprecationWarning('Please use atoms.calc')) + # def get_calculator(self): + # """Get currently attached calculator object. + + # Please use the equivalent atoms.calc instead of + # atoms.get_calculator().""" + # return self.calc + + # @property + # def calc(self): + # """Calculator object.""" + # return self._calc + + # @calc.setter + # def calc(self, calc): + # self._calc = calc + # if hasattr(calc, 'set_atoms'): + # calc.set_atoms(self) + + # @calc.deleter # type: ignore + # @deprecated(DeprecationWarning('Please use atoms.calc = None')) + # def calc(self): + # self._calc = None + + # @property # type: ignore + # @deprecated('Please use atoms.cell.rank instead') + # def number_of_lattice_vectors(self): + # """Number of (non-zero) lattice vectors.""" + # return self.cell.rank + + # def set_constraint(self, constraint=None): + # """Apply one or more constrains. + + # The *constraint* argument must be one constraint object or a + # list of constraint objects.""" + # if constraint is None: + # self._constraints = [] + # else: + # if isinstance(constraint, list): + # self._constraints = constraint + # elif isinstance(constraint, tuple): + # self._constraints = list(constraint) + # else: + # self._constraints = [constraint] + + # def _get_constraints(self): + # return self._constraints + + # def _del_constraints(self): + # self._constraints = [] + + # constraints = property( + # _get_constraints, set_constraint, _del_constraints, 'Constraints of the atoms.' + # ) def set_cell(self, cell, scale_atoms=False, apply_constraint=True): """Set unit cell vectors. @@ -381,7 +381,7 @@ def set_cell(self, cell, scale_atoms=False, apply_constraint=True): """ # Override pbcs if and only if given a Cell object: - cell = Cell.new(cell) + cell = ase.Cell.new(cell) # XXX not working well during initialize due to missing _constraints if apply_constraint and hasattr(self, '_constraints'): @@ -528,53 +528,6 @@ def has(self, name): # XXX extend has to calculator properties return name in self.arrays - def set_atomic_numbers(self, numbers): - """Set atomic numbers.""" - self.set_array('numbers', numbers, int, ()) - - def get_atomic_numbers(self): - """Get integer array of atomic numbers.""" - return self.arrays['numbers'].copy() - - def get_chemical_symbols(self): - """Get list of chemical symbol strings. - - Equivalent to ``list(atoms.symbols)``.""" - return list(self.symbols) - - def set_chemical_symbols(self, symbols): - """Set chemical symbols.""" - self.set_array('numbers', symbols2numbers(symbols), int, ()) - - def get_chemical_formula(self, mode='hill', empirical=False): - """Get the chemical formula as a string based on the chemical symbols. - - Parameters: - - mode: str - There are four different modes available: - - 'all': The list of chemical symbols are contracted to a string, - e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'. - - 'reduce': The same as 'all' where repeated elements are contracted - to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to - 'CH3OCH3'. - - 'hill': The list of chemical symbols are contracted to a string - following the Hill notation (alphabetical order with C and H - first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to - 'H2O4S'. This is default. - - 'metal': The list of chemical symbols (alphabetical metals, - and alphabetical non-metals) - - empirical, bool (optional, default=False) - Divide the symbol counts by their greatest common divisor to yield - an empirical formula. Only for mode `metal` and `hill`. - """ - return self.symbols.get_chemical_formula(mode, empirical) - def set_tags(self, tags): """Set tags for all atoms. If only one tag is supplied, it is applied to all atoms.""" @@ -609,79 +562,33 @@ def get_momenta(self): else: return np.zeros((len(self), 3)) - def set_masses(self, masses='defaults'): - """Set atomic masses in atomic mass units. - - The array masses should contain a list of masses. In case - the masses argument is not given or for those elements of the - masses list that are None, standard values are set.""" - - if isinstance(masses, str): - if masses == 'defaults': - masses = atomic_masses[self.arrays['numbers']] - elif masses == 'most_common': - masses = atomic_masses_common[self.arrays['numbers']] - elif masses is None: - pass - elif not isinstance(masses, np.ndarray): - masses = list(masses) - for i, mass in enumerate(masses): - if mass is None: - masses[i] = atomic_masses[self.numbers[i]] - self.set_array('masses', masses, float, ()) - - def get_masses(self): - """Get array of masses in atomic mass units.""" - if 'masses' in self.arrays: - return self.arrays['masses'].copy() - else: - return atomic_masses[self.arrays['numbers']] - - def set_initial_magnetic_moments(self, magmoms=None): - """Set the initial magnetic moments. - - Use either one or three numbers for every atom (collinear - or non-collinear spins).""" - - if magmoms is None: - self.set_array('initial_magmoms', None) - else: - magmoms = np.asarray(magmoms) - self.set_array('initial_magmoms', magmoms, float, magmoms.shape[1:]) - - def get_initial_magnetic_moments(self): - """Get array of initial magnetic moments.""" - if 'initial_magmoms' in self.arrays: - return self.arrays['initial_magmoms'].copy() - else: - return np.zeros(len(self)) - - def get_magnetic_moments(self): - """Get calculated local magnetic moments.""" - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - return self._calc.get_magnetic_moments(self) - - def get_magnetic_moment(self): - """Get calculated total magnetic moment.""" - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - return self._calc.get_magnetic_moment(self) - - def set_initial_charges(self, charges=None): - """Set the initial charges.""" - - if charges is None: - self.set_array('initial_charges', None) - else: - self.set_array('initial_charges', charges, float, ()) - - def get_initial_charges(self): - """Get array of initial charges.""" - if 'initial_charges' in self.arrays: - return self.arrays['initial_charges'].copy() - else: - return np.zeros(len(self)) + # def set_masses(self, masses='defaults'): + # """Set atomic masses in atomic mass units. + + # The array masses should contain a list of masses. In case + # the masses argument is not given or for those elements of the + # masses list that are None, standard values are set.""" + + # if isinstance(masses, str): + # if masses == 'defaults': + # masses = atomic_masses[self.arrays['numbers']] + # elif masses == 'most_common': + # masses = atomic_masses_common[self.arrays['numbers']] + # elif masses is None: + # pass + # elif not isinstance(masses, np.ndarray): + # masses = list(masses) + # for i, mass in enumerate(masses): + # if mass is None: + # masses[i] = atomic_masses[self.numbers[i]] + # self.set_array('masses', masses, float, ()) + + # def get_masses(self): + # """Get array of masses in atomic mass units.""" + # if 'masses' in self.arrays: + # return self.arrays['masses'].copy() + # else: + # return atomic_masses[self.arrays['numbers']] def get_charges(self): """Get calculated charges.""" @@ -704,48 +611,23 @@ def set_positions(self, newpositions, apply_constraint=True): self.set_array('positions', newpositions, shape=(3,)) - def get_positions(self, wrap=False, **wrap_kw): - """Get array of positions. + # def get_positions(self, wrap=False, **wrap_kw): + # """Get array of positions. - Parameters: + # Parameters: - wrap: bool - wrap atoms back to the cell before returning positions - wrap_kw: (keyword=value) pairs - optional keywords `pbc`, `center`, `pretty_translation`, `eps`, - see :func:`ase.geometry.wrap_positions` - """ - if wrap: - if 'pbc' not in wrap_kw: - wrap_kw['pbc'] = self.pbc - return wrap_positions(self.positions, self.cell, **wrap_kw) - else: - return self.arrays['positions'].copy() - - def get_potential_energy(self, force_consistent=False, apply_constraint=True): - """Calculate potential energy. - - Ask the attached calculator to calculate the potential energy and - apply constraints. Use *apply_constraint=False* to get the raw - forces. - - When supported by the calculator, either the energy extrapolated - to zero Kelvin or the energy consistent with the forces (the free - energy) can be returned. - """ - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - if force_consistent: - energy = self._calc.get_potential_energy( - self, force_consistent=force_consistent - ) - else: - energy = self._calc.get_potential_energy(self) - if apply_constraint: - for constraint in self.constraints: - if hasattr(constraint, 'adjust_potential_energy'): - energy += constraint.adjust_potential_energy(self) - return energy + # wrap: bool + # wrap atoms back to the cell before returning positions + # wrap_kw: (keyword=value) pairs + # optional keywords `pbc`, `center`, `pretty_translation`, `eps`, + # see :func:`ase.geometry.wrap_positions` + # """ + # if wrap: + # if 'pbc' not in wrap_kw: + # wrap_kw['pbc'] = self.pbc + # return wrap_positions(self.positions, self.cell, **wrap_kw) + # else: + # return self.arrays['positions'].copy() def get_properties(self, properties): """This method is experimental; currently for internal use.""" @@ -754,178 +636,23 @@ def get_properties(self, properties): raise RuntimeError('Atoms object has no calculator.') return self._calc.calculate_properties(self, properties) - def get_potential_energies(self): - """Calculate the potential energies of all the atoms. + # def get_velocities(self): + # """Get array of velocities.""" + # momenta = self.get_momenta() + # masses = self.get_masses() + # return momenta / masses[:, np.newaxis] - Only available with calculators supporting per-atom energies - (e.g. classical potentials). - """ - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - return self._calc.get_potential_energies(self) - - def get_kinetic_energy(self): - """Get the kinetic energy.""" - momenta = self.arrays.get('momenta') - if momenta is None: - return 0.0 - return 0.5 * np.vdot(momenta, self.get_velocities()) - - def get_velocities(self): - """Get array of velocities.""" - momenta = self.get_momenta() - masses = self.get_masses() - return momenta / masses[:, np.newaxis] - - def get_total_energy(self): - """Get the total energy - potential plus kinetic energy.""" - return self.get_potential_energy() + self.get_kinetic_energy() - - def get_forces(self, apply_constraint=True, md=False): - """Calculate atomic forces. - - Ask the attached calculator to calculate the forces and apply - constraints. Use *apply_constraint=False* to get the raw - forces. - - For molecular dynamics (md=True) we don't apply the constraint - to the forces but to the momenta. When holonomic constraints for - rigid linear triatomic molecules are present, ask the constraints - to redistribute the forces within each triple defined in the - constraints (required for molecular dynamics with this type of - constraints).""" - - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - forces = self._calc.get_forces(self) - - if apply_constraint: - # We need a special md flag here because for MD we want - # to skip real constraints but include special "constraints" - # Like Hookean. - for constraint in self.constraints: - if md and hasattr(constraint, 'redistribute_forces_md'): - constraint.redistribute_forces_md(self, forces) - if not md or hasattr(constraint, 'adjust_potential_energy'): - constraint.adjust_forces(self, forces) - return forces - - # Informs calculators (e.g. Asap) that ideal gas contribution is added here. - _ase_handles_dynamic_stress = True - - def get_stress(self, voigt=True, apply_constraint=True, include_ideal_gas=False): - """Calculate stress tensor. - - Returns an array of the six independent components of the - symmetric stress tensor, in the traditional Voigt order - (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt - order. - - The ideal gas contribution to the stresses is added if the - atoms have momenta and ``include_ideal_gas`` is set to True. - """ - - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - - stress = self._calc.get_stress(self) - shape = stress.shape - - if shape == (3, 3): - # Convert to the Voigt form before possibly applying - # constraints and adding the dynamic part of the stress - # (the "ideal gas contribution"). - stress = full_3x3_to_voigt_6_stress(stress) - else: - assert shape == (6,) - - if apply_constraint: - for constraint in self.constraints: - if hasattr(constraint, 'adjust_stress'): - constraint.adjust_stress(self, stress) - - # Add ideal gas contribution, if applicable - if include_ideal_gas and self.has('momenta'): - stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]]) - p = self.get_momenta() - masses = self.get_masses() - invmass = 1.0 / masses - invvol = 1.0 / self.get_volume() - for alpha in range(3): - for beta in range(alpha, 3): - stress[stresscomp[alpha, beta]] -= ( - p[:, alpha] * p[:, beta] * invmass - ).sum() * invvol - - if voigt: - return stress - else: - return voigt_6_to_full_3x3_stress(stress) - - def get_stresses(self, include_ideal_gas=False, voigt=True): - """Calculate the stress-tensor of all the atoms. - - Only available with calculators supporting per-atom energies and - stresses (e.g. classical potentials). Even for such calculators - there is a certain arbitrariness in defining per-atom stresses. - - The ideal gas contribution to the stresses is added if the - atoms have momenta and ``include_ideal_gas`` is set to True. - """ - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - stresses = self._calc.get_stresses(self) - - # make sure `stresses` are in voigt form - if np.shape(stresses)[1:] == (3, 3): - stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses] - stresses = np.array(stresses_voigt) - - # REMARK: The ideal gas contribution is intensive, i.e., the volume - # is divided out. We currently don't check if `stresses` are intensive - # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`. - # It might be good to check this here, but adds computational overhead. - - if include_ideal_gas and self.has('momenta'): - stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]]) - if hasattr(self._calc, 'get_atomic_volumes'): - invvol = 1.0 / self._calc.get_atomic_volumes() - else: - invvol = self.get_global_number_of_atoms() / self.get_volume() - p = self.get_momenta() - invmass = 1.0 / self.get_masses() - for alpha in range(3): - for beta in range(alpha, 3): - stresses[:, stresscomp[alpha, beta]] -= ( - p[:, alpha] * p[:, beta] * invmass * invvol - ) - if voigt: - return stresses - else: - stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses] - return np.array(stresses_3x3) - - def get_dipole_moment(self): - """Calculate the electric dipole moment for the atoms object. - - Only available for calculators which has a get_dipole_moment() - method.""" - - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - return self._calc.get_dipole_moment(self) - - def copy(self): - """Return a copy.""" - atoms = self.__class__( - cell=self.cell, pbc=self.pbc, info=self.info, celldisp=self._celldisp.copy() - ) + # def copy(self): + # """Return a copy.""" + # atoms = self.__class__( + # cell=self.cell, pbc=self.pbc, info=self.info, celldisp=self._celldisp.copy() + # ) - atoms.arrays = {} - for name, a in self.arrays.items(): - atoms.arrays[name] = a.copy() - atoms.constraints = copy.deepcopy(self.constraints) - return atoms + # atoms.arrays = {} + # for name, a in self.arrays.items(): + # atoms.arrays[name] = a.copy() + # atoms.constraints = copy.deepcopy(self.constraints) + # return atoms def todict(self): """For basic JSON (non-database) support.""" @@ -1047,7 +774,7 @@ def __add__(self, other): def extend(self, other): """Extend atoms object by appending atoms from *other*.""" - if isinstance(other, Atom): + if isinstance(other, Particle): other = self.__class__([other]) n1 = len(self) @@ -1107,7 +834,7 @@ def __getitem__(self, i): if i < -natoms or i >= natoms: raise IndexError('Index out of range.') - return Atom(atoms=self, index=i) + return Particle(atoms=self, index=i) elif not isinstance(i, slice): i = np.array(i) # if i is a mask @@ -1423,7 +1150,7 @@ def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False): a, v = v, a norm = np.linalg.norm - v = string2vector(v) + v = ase.string2vector(v) normv = norm(v) @@ -1431,12 +1158,12 @@ def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False): raise ZeroDivisionError('Cannot rotate: norm(v) == 0') if isinstance(a, numbers.Real): - a *= pi / 180 + a *= np.pi / 180 v /= normv - c = cos(a) - s = sin(a) + c = np.cos(a) + s = np.sin(a) else: - v2 = string2vector(a) + v2 = ase.string2vector(a) v /= normv normv2 = np.linalg.norm(v2) if normv2 == 0: @@ -1507,9 +1234,9 @@ def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)): """ center = self._centering_as_array(center) - phi *= pi / 180 - theta *= pi / 180 - psi *= pi / 180 + phi *= np.pi / 180 + theta *= np.pi / 180 + psi *= np.pi / 180 # First move the molecule to the origin In contrast to MATLAB, # numpy broadcasts the smaller array to the larger row-wise, @@ -1517,19 +1244,27 @@ def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)): rcoords = self.positions - center # First Euler rotation about z in matrix form D = np.array( - ((cos(phi), sin(phi), 0.0), (-sin(phi), cos(phi), 0.0), (0.0, 0.0, 1.0)) + ( + (np.cos(phi), np.sin(phi), 0.0), + (-np.sin(phi), np.cos(phi), 0.0), + (0.0, 0.0, 1.0), + ) ) # Second Euler rotation about x: C = np.array( ( (1.0, 0.0, 0.0), - (0.0, cos(theta), sin(theta)), - (0.0, -sin(theta), cos(theta)), + (0.0, np.cos(theta), np.sin(theta)), + (0.0, -np.sin(theta), np.cos(theta)), ) ) # Third Euler rotation, 2nd rotation about z: B = np.array( - ((cos(psi), sin(psi), 0.0), (-sin(psi), cos(psi), 0.0), (0.0, 0.0, 1.0)) + ( + (np.cos(psi), np.sin(psi), 0.0), + (-np.sin(psi), np.cos(psi), 0.0), + (0.0, 0.0, 1.0), + ) ) # Total Euler rotation A = np.dot(B, np.dot(C, D)) @@ -1578,7 +1313,7 @@ def get_dihedrals(self, indices, mic=False): cell = self.cell pbc = self.pbc - return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc) + return self.get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc) def _masked_rotate(self, center, axis, diff, mask): # do rotation of subgroup by copying it to temporary atoms object @@ -1591,7 +1326,7 @@ def _masked_rotate(self, center, axis, diff, mask): if mask[i]: group += self[i] group.translate(-center) - group.rotate(diff * 180 / pi, axis) + group.rotate(diff * 180 / np.pi, axis) group.translate(center) # set positions in original atoms object j = 0 @@ -1621,7 +1356,7 @@ def set_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None): >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1]) """ - angle *= pi / 180 + angle *= np.pi / 180 # if not provided, set mask to the last atom in the # dihedral description @@ -1632,7 +1367,7 @@ def set_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None): mask = [index in indices for index in range(len(self))] # compute necessary in dihedral change, from current value - current = self.get_dihedral(a1, a2, a3, a4) * pi / 180 + current = self.get_dihedral(a1, a2, a3, a4) * np.pi / 180 diff = angle - current axis = self.positions[a3] - self.positions[a2] center = self.positions[a3] @@ -1684,7 +1419,7 @@ def get_angles(self, indices, mic=False): cell = self.cell pbc = self.pbc - return get_angles(v12, v32, cell=cell, pbc=pbc) + return self.get_angles(v12, v32, cell=cell, pbc=pbc) def set_angle( self, a1, a2=None, a3=None, angle=None, mask=None, indices=None, add=False @@ -1716,7 +1451,7 @@ def set_angle( # Compute necessary in angle change, from current value diff = angle - self.get_angle(a1, a2, a3) - diff *= pi / 180 + diff *= np.pi / 180 # Do rotation of subgroup by copying it to temporary atoms object and # then rotating that v10 = self.positions[a1] - self.positions[a2] @@ -1727,26 +1462,6 @@ def set_angle( center = self.positions[a2] self._masked_rotate(center, axis, diff, mask) - def rattle(self, stdev=0.001, seed=None, rng=None): - """Randomly displace atoms. - - This method adds random displacements to the atomic positions, - taking a possible constraint into account. The random numbers are - drawn from a normal distribution of standard deviation stdev. - - For a parallel calculation, it is important to use the same - seed on all processors!""" - - if seed is not None and rng is not None: - raise ValueError('Please do not provide both seed and rng.') - - if rng is None: - if seed is None: - seed = 42 - rng = np.random.RandomState(seed) - positions = self.arrays['positions'] - self.set_positions(positions + rng.normal(scale=stdev, size=positions.shape)) - def get_distance(self, a0, a1, mic=False, vector=False): """Return distance between two atoms. @@ -1772,7 +1487,7 @@ def get_distances(self, a, indices, mic=False, vector=False): cell = self.cell pbc = self.pbc - D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc) + D, D_len = self.get_distances(p1, p2, cell=cell, pbc=pbc) if vector: D.shape = (-1, 3) @@ -1795,7 +1510,7 @@ def get_all_distances(self, mic=False, vector=False): cell = self.cell pbc = self.pbc - D, D_len = get_distances(R, cell=cell, pbc=pbc) + D, D_len = self.get_distances(R, cell=cell, pbc=pbc) if vector: return D @@ -1858,7 +1573,7 @@ def set_distance( D = np.array([R[a1] - R[a0]]) if mic: - D, D_len = find_mic(D, self.cell, self.pbc) + D, D_len = self.find_mic(D, self.cell, self.pbc) else: D_len = np.array([np.sqrt((D**2).sum())]) x = 1.0 - distance / D_len[0] @@ -1918,43 +1633,6 @@ def wrap(self, **wrap_kw): self.positions[:] = self.get_positions(wrap=True, **wrap_kw) - def get_temperature(self): - """Get the temperature in Kelvin.""" - dof = len(self) * 3 - for constraint in self._constraints: - dof -= constraint.get_removed_dof(self) - ekin = self.get_kinetic_energy() - return 2 * ekin / (dof * units.kB) - - def __eq__(self, other): - """Check for identity of two atoms objects. - - Identity means: same positions, atomic numbers, unit cell and - periodic boundary conditions.""" - if not isinstance(other, Atoms): - return False - a = self.arrays - b = other.arrays - return ( - len(self) == len(other) - and (a['positions'] == b['positions']).all() - and (a['numbers'] == b['numbers']).all() - and (self.cell == other.cell).all() - and (self.pbc == other.pbc).all() - ) - - def __ne__(self, other): - """Check if two atoms objects are not equal. - - Any differences in positions, atomic numbers, unit cell or - periodic boundary condtions make atoms objects not equal. - """ - eq = self.__eq__(other) - if eq is NotImplemented: - return eq - else: - return not eq - # @deprecated('Please use atoms.cell.volume') # We kind of want to deprecate this, but the ValueError behaviour # might be desirable. Should we do this? @@ -1980,87 +1658,6 @@ def _set_positions(self, pos): doc='Attribute for direct ' + 'manipulation of the positions.', ) - def _get_atomic_numbers(self): - """Return reference to atomic numbers for in-place - manipulations.""" - return self.arrays['numbers'] - - numbers = property( - _get_atomic_numbers, - set_atomic_numbers, - doc='Attribute for direct ' + 'manipulation of the atomic numbers.', - ) - - @property - def cell(self): - """The :class:`ase.cell.Cell` for direct manipulation.""" - return self._cellobj - - @cell.setter - def cell(self, cell): - cell = Cell.ascell(cell) - self._cellobj[:] = cell - - def write(self, filename, format=None, **kwargs): - """Write atoms object to a file. - - see ase.io.write for formats. - kwargs are passed to ase.io.write. - """ - from ase.io import write - - write(filename, self, format, **kwargs) - - def iterimages(self): - yield self - - def edit(self): - """Modify atoms interactively through ASE's GUI viewer. - - Conflicts leading to undesirable behaviour might arise - when matplotlib has been pre-imported with certain - incompatible backends and while trying to use the - plot feature inside the interactive GUI. To circumvent, - please set matplotlib.use('gtk') before calling this - method. - """ - from ase.gui.gui import GUI - from ase.gui.images import Images - - images = Images([self]) - gui = GUI(images) - gui.run() - - -def string2vector(v): - if isinstance(v, str): - if v[0] == '-': - return -string2vector(v[1:]) - w = np.zeros(3) - w['xyz'.index(v)] = 1.0 - return w - return np.array(v, float) - - -def default(data, dflt): - """Helper function for setting default values.""" - if data is None: - return None - elif isinstance(data, (list, tuple)): - newdata = [] - allnone = True - for x in data: - if x is None: - newdata.append(dflt) - else: - newdata.append(x) - allnone = False - if allnone: - return None - return newdata - else: - return data - # ? How generic (usable for any CG model) vs. Martini-specific do we want to be? class ParticlesState(Entity): From 9508fd975cbc5c5d0217095b214d624f7892e09d Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Mon, 14 Oct 2024 17:47:29 +0200 Subject: [PATCH 06/16] working on CG schema --- .../schema_packages/model_system.py | 188 +-- .../schema_packages/particles_state.py | 1405 ++--------------- 2 files changed, 192 insertions(+), 1401 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 4671bf56..75205d95 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -27,7 +27,7 @@ from structlog.stdlib import BoundLogger from nomad_simulations.schema_packages.atoms_state import AtomsState -from nomad_simulations.schema_packages.particles_state import ParticlesState +from nomad_simulations.schema_packages.particles_state import Particles, ParticlesState from nomad_simulations.schema_packages.utils import ( get_sibling_section, is_not_representative, @@ -553,99 +553,99 @@ def is_equal_cell(self, other) -> bool: return False return True - # def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]: - # """ - # Get the chemical symbols of the atoms in the atomic cell. These are defined on `atoms_state[*].chemical_symbol`. - - # Args: - # logger (BoundLogger): The logger to log messages. - - # Returns: - # list: The list of chemical symbols of the atoms in the atomic cell. - # """ - # if not self.atoms_state: - # return [] - - # chemical_symbols = [] - # for atom_state in self.atoms_state: - # if not atom_state.chemical_symbol: - # logger.warning('Could not find `AtomsState[*].chemical_symbol`.') - # return [] - # chemical_symbols.append(atom_state.chemical_symbol) - # return chemical_symbols - - # def to_ase_atoms(self, logger: 'BoundLogger') -> Optional[ase.Atoms]: - # """ - # Generates an ASE Atoms object with the most basic information from the parsed `AtomicCell` - # section (labels, periodic_boundary_conditions, positions, and lattice_vectors). - - # Args: - # logger (BoundLogger): The logger to log messages. - - # Returns: - # (Optional[ase.Atoms]): The ASE Atoms object with the basic information from the `AtomicCell`. - # """ - # # Initialize ase.Atoms object with labels - # atoms_labels = self.get_chemical_symbols(logger=logger) - # ase_atoms = ase.Atoms(symbols=atoms_labels) - - # # PBC - # if self.periodic_boundary_conditions is None: - # logger.info( - # 'Could not find `AtomicCell.periodic_boundary_conditions`. They will be set to [False, False, False].' - # ) - # self.periodic_boundary_conditions = [False, False, False] - # ase_atoms.set_pbc(pbc=self.periodic_boundary_conditions) - - # # Lattice vectors - # if self.lattice_vectors is not None: - # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) - # else: - # logger.info('Could not find `AtomicCell.lattice_vectors`.') - - # # Positions - # if self.positions is not None: - # if len(self.positions) != len(self.atoms_state): - # logger.error( - # 'Length of `AtomicCell.positions` does not coincide with the length of the `AtomicCell.atoms_state`.' - # ) - # return None - # ase_atoms.set_positions( - # newpositions=self.positions.to('angstrom').magnitude - # ) - # else: - # logger.warning('Could not find `AtomicCell.positions`.') - # return None - - # return ase_atoms - - # def from_ase_atoms(self, ase_atoms: ase.Atoms, logger: 'BoundLogger') -> None: - # """ - # Parses the information from an ASE Atoms object to the `AtomicCell` section. - - # Args: - # ase_atoms (ase.Atoms): The ASE Atoms object to parse. - # logger (BoundLogger): The logger to log messages. - # """ - # # `AtomsState[*].chemical_symbol` - # for symbol in ase_atoms.get_chemical_symbols(): - # atom_state = AtomsState(chemical_symbol=symbol) - # self.atoms_state.append(atom_state) - - # # `periodic_boundary_conditions` - # self.periodic_boundary_conditions = ase_atoms.get_pbc() - - # # `lattice_vectors` - # cell = ase_atoms.get_cell() - # self.lattice_vectors = ase.geometry.complete_cell(cell) * ureg('angstrom') - - # # `positions` - # positions = ase_atoms.get_positions() - # if ( - # not positions.tolist() - # ): # ASE assigns a shape=(0, 3) array if no positions are found - # return None - # self.positions = positions * ureg('angstrom') + def get_particle_types(self, logger: 'BoundLogger') -> list[str]: + """ + Get the chemical symbols of the atoms in the atomic cell. These are defined on `atoms_state[*].chemical_symbol`. + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + list: The list of chemical symbols of the atoms in the atomic cell. + """ + if not self.particles_state: + return [] + + particle_labels = [] + for particle_state in self.particles_state: + if not particle_state.particle_type: + logger.warning('Could not find `ParticlesState[*].particle_type`.') + return [] + particle_labels.append(particle_state.particle_type) + return particle_labels + + def to_particles(self, logger: 'BoundLogger') -> Optional[Particles]: + """ + Generates a Particles object with the most basic information from the parsed `ParticleCell` + section (labels, periodic_boundary_conditions, positions, and lattice_vectors). + + Args: + logger (BoundLogger): The logger to log messages. + + Returns: + (Optional[Particles]): The Partilces object with the basic information from the `ParticleCell`. + """ + # Initialize Partilces object with labels + particle_labels = self.get_particle_types(logger=logger) + particles = Particles(symbols=particle_labels) + + # PBC + if self.periodic_boundary_conditions is None: + logger.info( + 'Could not find `ParticleCell.periodic_boundary_conditions`. They will be set to [False, False, False].' + ) + self.periodic_boundary_conditions = [False, False, False] + particles.set_pbc(pbc=self.periodic_boundary_conditions) + + # Lattice vectors + if self.lattice_vectors is not None: + particles.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + else: + logger.info('Could not find `ParticleCell.lattice_vectors`.') + + # Positions + if self.positions is not None: + if len(self.positions) != len(self.particles_state): + logger.error( + 'Length of `ParticleCell.positions` does not coincide with the length of the `ParticleCell.particles_state`.' + ) + return None + particles.set_positions( + newpositions=self.positions.to('angstrom').magnitude + ) + else: + logger.warning('Could not find `ParticleCell.positions`.') + return None + + return particles + + def from_particles(self, particles: Particles, logger: 'BoundLogger') -> None: + """ + Parses the information from a Particles object to the `ParticlesCell` section. + + Args: + particles (Particles): The Particles object to parse. + logger (BoundLogger): The logger to log messages. + """ + # `ParticlesState[*].particles_type` + for label in particles.get_particle_types(): + particle_state = ParticlesState(particle_type=label) + self.particles_state.append(particle_state) + + # `periodic_boundary_conditions` + self.periodic_boundary_conditions = particles.get_pbc() + + # `lattice_vectors` + cell = particles.get_cell() + self.lattice_vectors = ase.geometry.complete_cell(cell) * ureg('angstrom') + + # `positions` + positions = particles.get_positions() + if ( + not positions.tolist() + ): # ASE assigns a shape=(0, 3) array if no positions are found + return None + self.positions = positions * ureg('angstrom') def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index b601e735..b4a9658a 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -2,6 +2,7 @@ from typing import TYPE_CHECKING, Any, Optional, Union import ase +import ase.geometry import numpy as np import pint from deprecated import deprecated @@ -17,20 +18,15 @@ from structlog.stdlib import BoundLogger -class Particle: +class Particles: """Particle object. - The Particle object can represent an isolated compound, or a - periodically repeated structure. It has a unit cell and - there may be periodic boundary conditions along any of the three - unit cell axes. - Information about the particles (types and position) is - stored in ndarrays. Optionally, there can be information about - type_shapes, mass, charge, diameter, body, moment_inertia, orientation, - angular momenta and image. - - In order to calculate energies, forces and stresses, a calculator - object has to attached to the atoms object. + Adaptation of the ASE Atoms object to coarse-grained particles. For use with + nomad_simulations.model_system.ParticlesCell. + Implemented methods: set_pbc, get_pbc, + set_cell, get_cell, + set_positions, get_positions, + get_particle_types Parameters: @@ -148,79 +144,79 @@ def __init__( particles = None - # if hasattr(symbols, 'get_positions'): - # atoms = symbols - # symbols = None - # elif ( - # isinstance(symbols, (list, tuple)) - # and len(symbols) > 0 - # and isinstance(symbols[0], Atom) - # ): - # # Get data from a list or tuple of Atom objects: - # data = [ - # [atom.get_raw(name) for atom in symbols] - # for name in [ - # 'position', - # 'number', - # 'tag', - # 'momentum', - # 'mass', - # 'magmom', - # 'charge', - # ] - # ] - # atoms = self.__class__(None, *data) - # symbols = None - - # if atoms is not None: - # # Get data from another Atoms object: - # if scaled_positions is not None: - # raise NotImplementedError - # if symbols is None and numbers is None: - # numbers = atoms.get_atomic_numbers() - # if positions is None: - # positions = atoms.get_positions() - # if tags is None and atoms.has('tags'): - # tags = atoms.get_tags() - # if momenta is None and atoms.has('momenta'): - # momenta = atoms.get_momenta() - # if magmoms is None and atoms.has('initial_magmoms'): - # magmoms = atoms.get_initial_magnetic_moments() - # if masses is None and atoms.has('masses'): - # masses = atoms.get_masses() - # if charges is None and atoms.has('initial_charges'): - # charges = atoms.get_initial_charges() - # if cell is None: - # cell = atoms.get_cell() - # if celldisp is None: - # celldisp = atoms.get_celldisp() - # if pbc is None: - # pbc = atoms.get_pbc() - - # self.arrays = {} - - # if symbols is None: - # if numbers is None: - # if positions is not None: - # natoms = len(positions) - # elif scaled_positions is not None: - # natoms = len(scaled_positions) - # else: - # natoms = 0 - # numbers = np.zeros(natoms, int) - # self.new_array('numbers', numbers, int) - # else: - # if numbers is not None: - # raise TypeError('Use only one of "symbols" and "numbers".') - # else: - # self.new_array('numbers', symbols2numbers(symbols), int) - - # if self.numbers.ndim != 1: - # raise ValueError('"numbers" must be 1-dimensional.') - - # if cell is None: - # cell = np.zeros((3, 3)) - # self.set_cell(cell) + # if hasattr(types, 'get_positions'): + # atoms = types + # types = None + # elif ( + # isinstance(types, (list, tuple)) + # and len(types) > 0 + # and isinstance(types[0], Atom) + # ): + # # Get data from a list or tuple of Atom objects: + # data = [ + # [atom.get_raw(name) for atom in types] + # for name in [ + # 'position', + # 'number', + # 'tag', + # 'momentum', + # 'mass', + # 'magmom', + # 'charge', + # ] + # ] + # atoms = self.__class__(None, *data) + # types = None + + # if atoms is not None: + # # Get data from another Atoms object: + # if scaled_positions is not None: + # raise NotImplementedError + # if types is None and numbers is None: + # numbers = atoms.get_atomic_numbers() + # if positions is None: + # positions = atoms.get_positions() + # if tags is None and atoms.has('tags'): + # tags = atoms.get_tags() + # if momenta is None and atoms.has('momenta'): + # momenta = atoms.get_momenta() + # if magmoms is None and atoms.has('initial_magmoms'): + # magmoms = atoms.get_initial_magnetic_moments() + # if masses is None and atoms.has('masses'): + # masses = atoms.get_masses() + # if charges is None and atoms.has('initial_charges'): + # charges = atoms.get_initial_charges() + # if cell is None: + # cell = atoms.get_cell() + # if celldisp is None: + # celldisp = atoms.get_celldisp() + # if pbc is None: + # pbc = atoms.get_pbc() + + # self.arrays = {} + + # if types is None: + # if numbers is None: + # if positions is not None: + # natoms = len(positions) + # elif scaled_positions is not None: + # natoms = len(scaled_positions) + # else: + # natoms = 0 + # numbers = np.zeros(natoms, int) + # self.new_array('numbers', numbers, int) + # else: + # if numbers is not None: + # raise TypeError('Use only one of "types" and "numbers".') + # else: + # self.new_array('numbers', types2numbers(types), int) + + # if self.numbers.ndim != 1: + # raise ValueError('"numbers" must be 1-dimensional.') + + if cell is None: + cell = np.zeros((3, 3)) + self.set_cell(cell) # if celldisp is None: # celldisp = np.zeros(shape=(3, 1)) @@ -234,7 +230,7 @@ def __init__( # positions = np.dot(scaled_positions, self.cell) # else: # if scaled_positions is not None: - # raise TypeError('Use only one of "symbols" and "numbers".') + # raise TypeError('Use only one of "types" and "numbers".') # self.new_array('positions', positions, float, (3,)) # self.set_tags(default(tags, 0)) # self.set_masses(default(masses, None)) @@ -275,68 +271,11 @@ def __init__( # new_symbols = Symbols.fromsymbols(obj) # self.numbers[:] = new_symbols.numbers - # @deprecated(DeprecationWarning('Please use atoms.calc = calc')) - # def set_calculator(self, calc=None): - # """Attach calculator object. - - # Please use the equivalent atoms.calc = calc instead of this - # method.""" - # self.calc = calc - - # @deprecated(DeprecationWarning('Please use atoms.calc')) - # def get_calculator(self): - # """Get currently attached calculator object. + def get_particle_types(self): + """Get list of particle type strings. - # Please use the equivalent atoms.calc instead of - # atoms.get_calculator().""" - # return self.calc - - # @property - # def calc(self): - # """Calculator object.""" - # return self._calc - - # @calc.setter - # def calc(self, calc): - # self._calc = calc - # if hasattr(calc, 'set_atoms'): - # calc.set_atoms(self) - - # @calc.deleter # type: ignore - # @deprecated(DeprecationWarning('Please use atoms.calc = None')) - # def calc(self): - # self._calc = None - - # @property # type: ignore - # @deprecated('Please use atoms.cell.rank instead') - # def number_of_lattice_vectors(self): - # """Number of (non-zero) lattice vectors.""" - # return self.cell.rank - - # def set_constraint(self, constraint=None): - # """Apply one or more constrains. - - # The *constraint* argument must be one constraint object or a - # list of constraint objects.""" - # if constraint is None: - # self._constraints = [] - # else: - # if isinstance(constraint, list): - # self._constraints = constraint - # elif isinstance(constraint, tuple): - # self._constraints = list(constraint) - # else: - # self._constraints = [constraint] - - # def _get_constraints(self): - # return self._constraints - - # def _del_constraints(self): - # self._constraints = [] - - # constraints = property( - # _get_constraints, set_constraint, _del_constraints, 'Constraints of the atoms.' - # ) + Labels describing type of coarse-grained particles.""" + return list(self.types) def set_cell(self, cell, scale_atoms=False, apply_constraint=True): """Set unit cell vectors. @@ -395,15 +334,6 @@ def set_cell(self, cell, scale_atoms=False, apply_constraint=True): self.cell[:] = cell - def set_celldisp(self, celldisp): - """Set the unit cell displacement vectors.""" - celldisp = np.array(celldisp, float) - self._celldisp = celldisp - - def get_celldisp(self): - """Get the unit cell displacement vectors.""" - return self._celldisp.copy() - def get_cell(self, complete=False): """Get the three unit cell vectors as a `class`:ase.cell.Cell` object. @@ -416,28 +346,6 @@ def get_cell(self, complete=False): return cell - @deprecated('Please use atoms.cell.cellpar() instead') - def get_cell_lengths_and_angles(self): - """Get unit cell parameters. Sequence of 6 numbers. - - First three are unit cell vector lengths and second three - are angles between them:: - - [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)] - - in degrees. - """ - return self.cell.cellpar() - - @deprecated('Please use atoms.cell.reciprocal()') - def get_reciprocal_cell(self): - """Get the three reciprocal lattice vectors as a 3x3 ndarray. - - Note that the commonly used factor of 2 pi for Fourier - transforms is not included here.""" - - return self.cell.reciprocal() - @property def pbc(self): """Reference to pbc-flags for in-place manipulations.""" @@ -455,152 +363,6 @@ def get_pbc(self): """Get periodic boundary condition flags.""" return self.pbc.copy() - def new_array(self, name, a, dtype=None, shape=None): - """Add new array. - - If *shape* is not *None*, the shape of *a* will be checked.""" - - if dtype is not None: - a = np.array(a, dtype, order='C') - if len(a) == 0 and shape is not None: - a.shape = (-1,) + shape - else: - if not a.flags['C_CONTIGUOUS']: - a = np.ascontiguousarray(a) - else: - a = a.copy() - - if name in self.arrays: - raise RuntimeError(f'Array {name} already present') - - for b in self.arrays.values(): - if len(a) != len(b): - raise ValueError( - 'Array "%s" has wrong length: %d != %d.' % (name, len(a), len(b)) - ) - break - - if shape is not None and a.shape[1:] != shape: - raise ValueError( - 'Array "%s" has wrong shape %s != %s.' - % (name, a.shape, (a.shape[0:1] + shape)) - ) - - self.arrays[name] = a - - def get_array(self, name, copy=True): - """Get an array. - - Returns a copy unless the optional argument copy is false. - """ - if copy: - return self.arrays[name].copy() - else: - return self.arrays[name] - - def set_array(self, name, a, dtype=None, shape=None): - """Update array. - - If *shape* is not *None*, the shape of *a* will be checked. - If *a* is *None*, then the array is deleted.""" - - b = self.arrays.get(name) - if b is None: - if a is not None: - self.new_array(name, a, dtype, shape) - else: - if a is None: - del self.arrays[name] - else: - a = np.asarray(a) - if a.shape != b.shape: - raise ValueError( - 'Array "%s" has wrong shape %s != %s.' - % (name, a.shape, b.shape) - ) - b[:] = a - - def has(self, name): - """Check for existence of array. - - name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms', - 'initial_charges'.""" - # XXX extend has to calculator properties - return name in self.arrays - - def set_tags(self, tags): - """Set tags for all atoms. If only one tag is supplied, it is - applied to all atoms.""" - if isinstance(tags, int): - tags = [tags] * len(self) - self.set_array('tags', tags, int, ()) - - def get_tags(self): - """Get integer array of tags.""" - if 'tags' in self.arrays: - return self.arrays['tags'].copy() - else: - return np.zeros(len(self), int) - - def set_momenta(self, momenta, apply_constraint=True): - """Set momenta.""" - if apply_constraint and len(self.constraints) > 0 and momenta is not None: - momenta = np.array(momenta) # modify a copy - for constraint in self.constraints: - if hasattr(constraint, 'adjust_momenta'): - constraint.adjust_momenta(self, momenta) - self.set_array('momenta', momenta, float, (3,)) - - def set_velocities(self, velocities): - """Set the momenta by specifying the velocities.""" - self.set_momenta(self.get_masses()[:, np.newaxis] * velocities) - - def get_momenta(self): - """Get array of momenta.""" - if 'momenta' in self.arrays: - return self.arrays['momenta'].copy() - else: - return np.zeros((len(self), 3)) - - # def set_masses(self, masses='defaults'): - # """Set atomic masses in atomic mass units. - - # The array masses should contain a list of masses. In case - # the masses argument is not given or for those elements of the - # masses list that are None, standard values are set.""" - - # if isinstance(masses, str): - # if masses == 'defaults': - # masses = atomic_masses[self.arrays['numbers']] - # elif masses == 'most_common': - # masses = atomic_masses_common[self.arrays['numbers']] - # elif masses is None: - # pass - # elif not isinstance(masses, np.ndarray): - # masses = list(masses) - # for i, mass in enumerate(masses): - # if mass is None: - # masses[i] = atomic_masses[self.numbers[i]] - # self.set_array('masses', masses, float, ()) - - # def get_masses(self): - # """Get array of masses in atomic mass units.""" - # if 'masses' in self.arrays: - # return self.arrays['masses'].copy() - # else: - # return atomic_masses[self.arrays['numbers']] - - def get_charges(self): - """Get calculated charges.""" - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - try: - return self._calc.get_charges(self) - except AttributeError: - from ase.calculators.calculator import PropertyNotImplementedError - - raise PropertyNotImplementedError - def set_positions(self, newpositions, apply_constraint=True): """Set positions, honoring any constraints. To ignore constraints, use *apply_constraint=False*.""" @@ -611,983 +373,23 @@ def set_positions(self, newpositions, apply_constraint=True): self.set_array('positions', newpositions, shape=(3,)) - # def get_positions(self, wrap=False, **wrap_kw): - # """Get array of positions. - - # Parameters: - - # wrap: bool - # wrap atoms back to the cell before returning positions - # wrap_kw: (keyword=value) pairs - # optional keywords `pbc`, `center`, `pretty_translation`, `eps`, - # see :func:`ase.geometry.wrap_positions` - # """ - # if wrap: - # if 'pbc' not in wrap_kw: - # wrap_kw['pbc'] = self.pbc - # return wrap_positions(self.positions, self.cell, **wrap_kw) - # else: - # return self.arrays['positions'].copy() - - def get_properties(self, properties): - """This method is experimental; currently for internal use.""" - # XXX Something about constraints. - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - return self._calc.calculate_properties(self, properties) - - # def get_velocities(self): - # """Get array of velocities.""" - # momenta = self.get_momenta() - # masses = self.get_masses() - # return momenta / masses[:, np.newaxis] - - # def copy(self): - # """Return a copy.""" - # atoms = self.__class__( - # cell=self.cell, pbc=self.pbc, info=self.info, celldisp=self._celldisp.copy() - # ) - - # atoms.arrays = {} - # for name, a in self.arrays.items(): - # atoms.arrays[name] = a.copy() - # atoms.constraints = copy.deepcopy(self.constraints) - # return atoms - - def todict(self): - """For basic JSON (non-database) support.""" - d = dict(self.arrays) - d['cell'] = np.asarray(self.cell) - d['pbc'] = self.pbc - if self._celldisp.any(): - d['celldisp'] = self._celldisp - if self.constraints: - d['constraints'] = self.constraints - if self.info: - d['info'] = self.info - # Calculator... trouble. - return d - - @classmethod - def fromdict(cls, dct): - """Rebuild atoms object from dictionary representation (todict).""" - dct = dct.copy() - kw = {} - for name in ['numbers', 'positions', 'cell', 'pbc']: - kw[name] = dct.pop(name) - - constraints = dct.pop('constraints', None) - if constraints: - from ase.constraints import dict2constraint - - constraints = [dict2constraint(d) for d in constraints] - - info = dct.pop('info', None) - - atoms = cls( - constraint=constraints, celldisp=dct.pop('celldisp', None), info=info, **kw - ) - natoms = len(atoms) - - # Some arrays are named differently from the atoms __init__ keywords. - # Also, there may be custom arrays. Hence we set them directly: - for name, arr in dct.items(): - assert len(arr) == natoms, name - assert isinstance(arr, np.ndarray) - atoms.arrays[name] = arr - return atoms - - def __len__(self): - return len(self.arrays['positions']) - - def get_number_of_atoms(self): - """Deprecated, please do not use. - - You probably want len(atoms). Or if your atoms are distributed, - use (and see) get_global_number_of_atoms().""" - import warnings - - warnings.warn( - 'Use get_global_number_of_atoms() instead', np.VisibleDeprecationWarning - ) - return len(self) - - def get_global_number_of_atoms(self): - """Returns the global number of atoms in a distributed-atoms parallel - simulation. - - DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING! - - Equivalent to len(atoms) in the standard ASE Atoms class. You should - normally use len(atoms) instead. This function's only purpose is to - make compatibility between ASE and Asap easier to maintain by having a - few places in ASE use this function instead. It is typically only - when counting the global number of degrees of freedom or in similar - situations. - """ - return len(self) - - def __repr__(self): - tokens = [] - - N = len(self) - if N <= 60: - symbols = self.get_chemical_formula('reduce') - else: - symbols = self.get_chemical_formula('hill') - tokens.append(f"symbols='{symbols}'") - - if self.pbc.any() and not self.pbc.all(): - tokens.append(f'pbc={self.pbc.tolist()}') - else: - tokens.append(f'pbc={self.pbc[0]}') - - cell = self.cell - if cell: - if cell.orthorhombic: - cell = cell.lengths().tolist() - else: - cell = cell.tolist() - tokens.append(f'cell={cell}') - - for name in sorted(self.arrays): - if name in ['numbers', 'positions']: - continue - tokens.append(f'{name}=...') - - if self.constraints: - if len(self.constraints) == 1: - constraint = self.constraints[0] - else: - constraint = self.constraints - tokens.append(f'constraint={repr(constraint)}') - - if self._calc is not None: - tokens.append(f'calculator={self._calc.__class__.__name__}(...)') - - return '{0}({1})'.format(self.__class__.__name__, ', '.join(tokens)) - - def __add__(self, other): - atoms = self.copy() - atoms += other - return atoms - - def extend(self, other): - """Extend atoms object by appending atoms from *other*.""" - if isinstance(other, Particle): - other = self.__class__([other]) - - n1 = len(self) - n2 = len(other) - - for name, a1 in self.arrays.items(): - a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype) - a[:n1] = a1 - if name == 'masses': - a2 = other.get_masses() - else: - a2 = other.arrays.get(name) - if a2 is not None: - a[n1:] = a2 - self.arrays[name] = a - - for name, a2 in other.arrays.items(): - if name in self.arrays: - continue - a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype) - a[n1:] = a2 - if name == 'masses': - a[:n1] = self.get_masses()[:n1] - else: - a[:n1] = 0 - - self.set_array(name, a) - - def __iadd__(self, other): - self.extend(other) - return self - - def append(self, atom): - """Append atom to end.""" - self.extend(self.__class__([atom])) - - def __iter__(self): - for i in range(len(self)): - yield self[i] - - def __getitem__(self, i): - """Return a subset of the atoms. - - i -- scalar integer, list of integers, or slice object - describing which atoms to return. - - If i is a scalar, return an Atom object. If i is a list or a - slice, return an Atoms object with the same cell, pbc, and - other associated info as the original Atoms object. The - indices of the constraints will be shuffled so that they match - the indexing in the subset returned. - - """ - - if isinstance(i, numbers.Integral): - natoms = len(self) - if i < -natoms or i >= natoms: - raise IndexError('Index out of range.') - - return Particle(atoms=self, index=i) - elif not isinstance(i, slice): - i = np.array(i) - # if i is a mask - if i.dtype == bool: - if len(i) != len(self): - raise IndexError( - f'Length of mask {len(i)} must equal ' - f'number of atoms {len(self)}' - ) - i = np.arange(len(self))[i] - - import copy - - conadd = [] - # Constraints need to be deepcopied, but only the relevant ones. - for con in copy.deepcopy(self.constraints): - try: - con.index_shuffle(self, i) - except (IndexError, NotImplementedError): - pass - else: - conadd.append(con) - - atoms = self.__class__( - cell=self.cell, - pbc=self.pbc, - info=self.info, - # should be communicated to the slice as well - celldisp=self._celldisp, - ) - # TODO: Do we need to shuffle indices in adsorbate_info too? - - atoms.arrays = {} - for name, a in self.arrays.items(): - atoms.arrays[name] = a[i].copy() - - atoms.constraints = conadd - return atoms - - def __delitem__(self, i): - from ase.constraints import FixAtoms - - for c in self._constraints: - if not isinstance(c, FixAtoms): - raise RuntimeError( - 'Remove constraint using set_constraint() ' 'before deleting atoms.' - ) - - if isinstance(i, list) and len(i) > 0: - # Make sure a list of booleans will work correctly and not be - # interpreted at 0 and 1 indices. - i = np.array(i) - - if len(self._constraints) > 0: - n = len(self) - i = np.arange(n)[i] - if isinstance(i, int): - i = [i] - constraints = [] - for c in self._constraints: - c = c.delete_atoms(i, n) - if c is not None: - constraints.append(c) - self.constraints = constraints - - mask = np.ones(len(self), bool) - mask[i] = False - for name, a in self.arrays.items(): - self.arrays[name] = a[mask] - - def pop(self, i=-1): - """Remove and return atom at index *i* (default last).""" - atom = self[i] - atom.cut_reference_to_atoms() - del self[i] - return atom - - def __imul__(self, m): - """In-place repeat of atoms.""" - if isinstance(m, int): - m = (m, m, m) - - for x, vec in zip(m, self.cell): - if x != 1 and not vec.any(): - raise ValueError('Cannot repeat along undefined lattice ' 'vector') - - M = np.product(m) - n = len(self) - - for name, a in self.arrays.items(): - self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1)) - - positions = self.arrays['positions'] - i0 = 0 - for m0 in range(m[0]): - for m1 in range(m[1]): - for m2 in range(m[2]): - i1 = i0 + n - positions[i0:i1] += np.dot((m0, m1, m2), self.cell) - i0 = i1 - - if self.constraints is not None: - self.constraints = [c.repeat(m, n) for c in self.constraints] - - self.cell = np.array([m[c] * self.cell[c] for c in range(3)]) - - return self - - def repeat(self, rep): - """Create new repeated atoms object. - - The *rep* argument should be a sequence of three positive - integers like *(2,3,1)* or a single integer (*r*) equivalent - to *(r,r,r)*.""" - - atoms = self.copy() - atoms *= rep - return atoms - - def __mul__(self, rep): - return self.repeat(rep) - - def translate(self, displacement): - """Translate atomic positions. - - The displacement argument can be a float an xyz vector or an - nx3 array (where n is the number of atoms).""" - - self.arrays['positions'] += np.array(displacement) - - def center(self, vacuum=None, axis=(0, 1, 2), about=None): - """Center atoms in unit cell. - - Centers the atoms in the unit cell, so there is the same - amount of vacuum on all sides. - - vacuum: float (default: None) - If specified adjust the amount of vacuum when centering. - If vacuum=10.0 there will thus be 10 Angstrom of vacuum - on each side. - axis: int or sequence of ints - Axis or axes to act on. Default: Act on all axes. - about: float or array (default: None) - If specified, center the atoms about . - I.e., about=(0., 0., 0.) (or just "about=0.", interpreted - identically), to center about the origin. - """ - - # Find the orientations of the faces of the unit cell - cell = self.cell.complete() - dirs = np.zeros_like(cell) - - lengths = cell.lengths() - for i in range(3): - dirs[i] = np.cross(cell[i - 1], cell[i - 2]) - dirs[i] /= np.linalg.norm(dirs[i]) - if dirs[i] @ cell[i] < 0.0: - dirs[i] *= -1 - - if isinstance(axis, int): - axes = (axis,) - else: - axes = axis - - # Now, decide how much each basis vector should be made longer - pos = self.positions - longer = np.zeros(3) - shift = np.zeros(3) - for i in axes: - if len(pos): - scalarprod = pos @ dirs[i] - p0 = scalarprod.min() - p1 = scalarprod.max() - else: - p0 = 0 - p1 = 0 - height = cell[i] @ dirs[i] - if vacuum is not None: - lng = (p1 - p0 + 2 * vacuum) - height - else: - lng = 0.0 # Do not change unit cell size! - top = lng + height - p1 - shf = 0.5 * (top - p0) - cosphi = cell[i] @ dirs[i] / lengths[i] - longer[i] = lng / cosphi - shift[i] = shf / cosphi - - # Now, do it! - translation = np.zeros(3) - for i in axes: - nowlen = lengths[i] - if vacuum is not None: - self.cell[i] = cell[i] * (1 + longer[i] / nowlen) - translation += shift[i] * cell[i] / nowlen - - # We calculated translations using the completed cell, - # so directions without cell vectors will have been centered - # along a "fake" vector of length 1. - # Therefore, we adjust by -0.5: - if not any(self.cell[i]): - translation[i] -= 0.5 - - # Optionally, translate to center about a point in space. - if about is not None: - for vector in self.cell: - translation -= vector / 2.0 - translation += about - - self.positions += translation - - def get_center_of_mass(self, scaled=False): - """Get the center of mass. - - If scaled=True the center of mass in scaled coordinates - is returned.""" - masses = self.get_masses() - com = masses @ self.positions / masses.sum() - if scaled: - return self.cell.scaled_positions(com) - else: - return com - - def set_center_of_mass(self, com, scaled=False): - """Set the center of mass. - - If scaled=True the center of mass is expected in scaled coordinates. - Constraints are considered for scaled=False. - """ - old_com = self.get_center_of_mass(scaled=scaled) - difference = old_com - com - if scaled: - self.set_scaled_positions(self.get_scaled_positions() + difference) - else: - self.set_positions(self.get_positions() + difference) - - def get_moments_of_inertia(self, vectors=False): - """Get the moments of inertia along the principal axes. - - The three principal moments of inertia are computed from the - eigenvalues of the symmetric inertial tensor. Periodic boundary - conditions are ignored. Units of the moments of inertia are - amu*angstrom**2. - """ - com = self.get_center_of_mass() - positions = self.get_positions() - positions -= com # translate center of mass to origin - masses = self.get_masses() - - # Initialize elements of the inertial tensor - I11 = I22 = I33 = I12 = I13 = I23 = 0.0 - for i in range(len(self)): - x, y, z = positions[i] - m = masses[i] - - I11 += m * (y**2 + z**2) - I22 += m * (x**2 + z**2) - I33 += m * (x**2 + y**2) - I12 += -m * x * y - I13 += -m * x * z - I23 += -m * y * z - - I = np.array([[I11, I12, I13], [I12, I22, I23], [I13, I23, I33]]) - - evals, evecs = np.linalg.eigh(I) - if vectors: - return evals, evecs.transpose() - else: - return evals - - def get_angular_momentum(self): - """Get total angular momentum with respect to the center of mass.""" - com = self.get_center_of_mass() - positions = self.get_positions() - positions -= com # translate center of mass to origin - return np.cross(positions, self.get_momenta()).sum(0) - - def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False): - """Rotate atoms based on a vector and an angle, or two vectors. - - Parameters: - - a = None: - Angle that the atoms is rotated around the vector 'v'. 'a' - can also be a vector and then 'a' is rotated - into 'v'. - - v: - Vector to rotate the atoms around. Vectors can be given as - strings: 'x', '-x', 'y', ... . - - center = (0, 0, 0): - The center is kept fixed under the rotation. Use 'COM' to fix - the center of mass, 'COP' to fix the center of positions or - 'COU' to fix the center of cell. - - rotate_cell = False: - If true the cell is also rotated. - - Examples: - - Rotate 90 degrees around the z-axis, so that the x-axis is - rotated into the y-axis: - - >>> atoms = Atoms() - >>> atoms.rotate(90, 'z') - >>> atoms.rotate(90, (0, 0, 1)) - >>> atoms.rotate(-90, '-z') - >>> atoms.rotate('x', 'y') - >>> atoms.rotate((1, 0, 0), (0, 1, 0)) - """ - - if not isinstance(a, numbers.Real): - a, v = v, a - - norm = np.linalg.norm - v = ase.string2vector(v) - - normv = norm(v) - - if normv == 0.0: - raise ZeroDivisionError('Cannot rotate: norm(v) == 0') - - if isinstance(a, numbers.Real): - a *= np.pi / 180 - v /= normv - c = np.cos(a) - s = np.sin(a) - else: - v2 = ase.string2vector(a) - v /= normv - normv2 = np.linalg.norm(v2) - if normv2 == 0: - raise ZeroDivisionError('Cannot rotate: norm(a) == 0') - v2 /= norm(v2) - c = np.dot(v, v2) - v = np.cross(v, v2) - s = norm(v) - # In case *v* and *a* are parallel, np.cross(v, v2) vanish - # and can't be used as a rotation axis. However, in this - # case any rotation axis perpendicular to v2 will do. - eps = 1e-7 - if s < eps: - v = np.cross((0, 0, 1), v2) - if norm(v) < eps: - v = np.cross((1, 0, 0), v2) - assert norm(v) >= eps - elif s > 0: - v /= s - - center = self._centering_as_array(center) - - p = self.arrays['positions'] - center - self.arrays['positions'][:] = ( - c * p - np.cross(p, s * v) + np.outer(np.dot(p, v), (1.0 - c) * v) + center - ) - if rotate_cell: - rotcell = self.get_cell() - rotcell[:] = ( - c * rotcell - - np.cross(rotcell, s * v) - + np.outer(np.dot(rotcell, v), (1.0 - c) * v) - ) - self.set_cell(rotcell) - - def _centering_as_array(self, center): - if isinstance(center, str): - if center.lower() == 'com': - center = self.get_center_of_mass() - elif center.lower() == 'cop': - center = self.get_positions().mean(axis=0) - elif center.lower() == 'cou': - center = self.get_cell().sum(axis=0) / 2 - else: - raise ValueError('Cannot interpret center') - else: - center = np.array(center, float) - return center - - def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)): - """Rotate atoms via Euler angles (in degrees). - - See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation. + def get_positions(self, wrap=False, **wrap_kw): + """Get array of positions. Parameters: - center : - The point to rotate about. A sequence of length 3 with the - coordinates, or 'COM' to select the center of mass, 'COP' to - select center of positions or 'COU' to select center of cell. - phi : - The 1st rotation angle around the z axis. - theta : - Rotation around the x axis. - psi : - 2nd rotation around the z axis. - - """ - center = self._centering_as_array(center) - - phi *= np.pi / 180 - theta *= np.pi / 180 - psi *= np.pi / 180 - - # First move the molecule to the origin In contrast to MATLAB, - # numpy broadcasts the smaller array to the larger row-wise, - # so there is no need to play with the Kronecker product. - rcoords = self.positions - center - # First Euler rotation about z in matrix form - D = np.array( - ( - (np.cos(phi), np.sin(phi), 0.0), - (-np.sin(phi), np.cos(phi), 0.0), - (0.0, 0.0, 1.0), - ) - ) - # Second Euler rotation about x: - C = np.array( - ( - (1.0, 0.0, 0.0), - (0.0, np.cos(theta), np.sin(theta)), - (0.0, -np.sin(theta), np.cos(theta)), - ) - ) - # Third Euler rotation, 2nd rotation about z: - B = np.array( - ( - (np.cos(psi), np.sin(psi), 0.0), - (-np.sin(psi), np.cos(psi), 0.0), - (0.0, 0.0, 1.0), - ) - ) - # Total Euler rotation - A = np.dot(B, np.dot(C, D)) - # Do the rotation - rcoords = np.dot(A, np.transpose(rcoords)) - # Move back to the rotation point - self.positions = np.transpose(rcoords) + center - - def get_dihedral(self, a0, a1, a2, a3, mic=False): - """Calculate dihedral angle. - - Calculate dihedral angle (in degrees) between the vectors a0->a1 - and a2->a3. - - Use mic=True to use the Minimum Image Convention and calculate the - angle across periodic boundaries. - """ - return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0] - - def get_dihedrals(self, indices, mic=False): - """Calculate dihedral angles. - - Calculate dihedral angles (in degrees) between the list of vectors - a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices. - - Use mic=True to use the Minimum Image Convention and calculate the - angles across periodic boundaries. - """ - indices = np.array(indices) - assert indices.shape[1] == 4 - - a0s = self.positions[indices[:, 0]] - a1s = self.positions[indices[:, 1]] - a2s = self.positions[indices[:, 2]] - a3s = self.positions[indices[:, 3]] - - # vectors 0->1, 1->2, 2->3 - v0 = a1s - a0s - v1 = a2s - a1s - v2 = a3s - a2s - - cell = None - pbc = None - - if mic: - cell = self.cell - pbc = self.pbc - - return self.get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc) - - def _masked_rotate(self, center, axis, diff, mask): - # do rotation of subgroup by copying it to temporary atoms object - # and then rotating that - # - # recursive object definition might not be the most elegant thing, - # more generally useful might be a rotation function with a mask? - group = self.__class__() - for i in range(len(self)): - if mask[i]: - group += self[i] - group.translate(-center) - group.rotate(diff * 180 / np.pi, axis) - group.translate(center) - # set positions in original atoms object - j = 0 - for i in range(len(self)): - if mask[i]: - self.positions[i] = group[j].position - j += 1 - - def set_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None): - """Set the dihedral angle (degrees) between vectors a1->a2 and - a3->a4 by changing the atom indexed by a4. - - If mask is not None, all the atoms described in mask - (read: the entire subgroup) are moved. Alternatively to the mask, - the indices of the atoms to be rotated can be supplied. If both - *mask* and *indices* are given, *indices* overwrites *mask*. - - **Important**: If *mask* or *indices* is given and does not contain - *a4*, *a4* will NOT be moved. In most cases you therefore want - to include *a4* in *mask*/*indices*. - - Example: the following defines a very crude - ethane-like molecule and twists one half of it by 30 degrees. - - >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0], - ... [1, 0, 0], [2, 1, 0], [2, -1, 0]]) - >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1]) - """ - - angle *= np.pi / 180 - - # if not provided, set mask to the last atom in the - # dihedral description - if mask is None and indices is None: - mask = np.zeros(len(self)) - mask[a4] = 1 - elif indices is not None: - mask = [index in indices for index in range(len(self))] - - # compute necessary in dihedral change, from current value - current = self.get_dihedral(a1, a2, a3, a4) * np.pi / 180 - diff = angle - current - axis = self.positions[a3] - self.positions[a2] - center = self.positions[a3] - self._masked_rotate(center, axis, diff, mask) - - def rotate_dihedral(self, a1, a2, a3, a4, angle=None, mask=None, indices=None): - """Rotate dihedral angle. - - Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a - predefined dihedral angle, starting from its current configuration. - """ - start = self.get_dihedral(a1, a2, a3, a4) - self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices) - - def get_angle(self, a1, a2, a3, mic=False): - """Get angle formed by three atoms. - - Calculate angle in degrees between the vectors a2->a1 and - a2->a3. - - Use mic=True to use the Minimum Image Convention and calculate the - angle across periodic boundaries. - """ - return self.get_angles([[a1, a2, a3]], mic=mic)[0] - - def get_angles(self, indices, mic=False): - """Get angle formed by three atoms for multiple groupings. - - Calculate angle in degrees between vectors between atoms a2->a1 - and a2->a3, where a1, a2, and a3 are in each row of indices. - - Use mic=True to use the Minimum Image Convention and calculate - the angle across periodic boundaries. - """ - indices = np.array(indices) - assert indices.shape[1] == 3 - - a1s = self.positions[indices[:, 0]] - a2s = self.positions[indices[:, 1]] - a3s = self.positions[indices[:, 2]] - - v12 = a1s - a2s - v32 = a3s - a2s - - cell = None - pbc = None - - if mic: - cell = self.cell - pbc = self.pbc - - return self.get_angles(v12, v32, cell=cell, pbc=pbc) - - def set_angle( - self, a1, a2=None, a3=None, angle=None, mask=None, indices=None, add=False - ): - """Set angle (in degrees) formed by three atoms. - - Sets the angle between vectors *a2*->*a1* and *a2*->*a3*. - - If *add* is `True`, the angle will be changed by the value given. - - Same usage as in :meth:`ase.Atoms.set_dihedral`. - If *mask* and *indices* - are given, *indices* overwrites *mask*. If *mask* and *indices* - are not set, only *a3* is moved.""" - - if any(a is None for a in [a2, a3, angle]): - raise ValueError('a2, a3, and angle must not be None') - - # If not provided, set mask to the last atom in the angle description - if mask is None and indices is None: - mask = np.zeros(len(self)) - mask[a3] = 1 - elif indices is not None: - mask = [index in indices for index in range(len(self))] - - if add: - diff = angle - else: - # Compute necessary in angle change, from current value - diff = angle - self.get_angle(a1, a2, a3) - - diff *= np.pi / 180 - # Do rotation of subgroup by copying it to temporary atoms object and - # then rotating that - v10 = self.positions[a1] - self.positions[a2] - v12 = self.positions[a3] - self.positions[a2] - v10 /= np.linalg.norm(v10) - v12 /= np.linalg.norm(v12) - axis = np.cross(v10, v12) - center = self.positions[a2] - self._masked_rotate(center, axis, diff, mask) - - def get_distance(self, a0, a1, mic=False, vector=False): - """Return distance between two atoms. - - Use mic=True to use the Minimum Image Convention. - vector=True gives the distance vector (from a0 to a1). - """ - return self.get_distances(a0, [a1], mic=mic, vector=vector)[0] - - def get_distances(self, a, indices, mic=False, vector=False): - """Return distances of atom No.i with a list of atoms. - - Use mic=True to use the Minimum Image Convention. - vector=True gives the distance vector (from a to self[indices]). - """ - R = self.arrays['positions'] - p1 = [R[a]] - p2 = R[indices] - - cell = None - pbc = None - - if mic: - cell = self.cell - pbc = self.pbc - - D, D_len = self.get_distances(p1, p2, cell=cell, pbc=pbc) - - if vector: - D.shape = (-1, 3) - return D - else: - D_len.shape = (-1,) - return D_len - - def get_all_distances(self, mic=False, vector=False): - """Return distances of all of the atoms with all of the atoms. - - Use mic=True to use the Minimum Image Convention. + wrap: bool + wrap atoms back to the cell before returning positions + wrap_kw: (keyword=value) pairs + optional keywords `pbc`, `center`, `pretty_translation`, `eps`, + see :func:`ase.geometry.wrap_positions` """ - R = self.arrays['positions'] - - cell = None - pbc = None - - if mic: - cell = self.cell - pbc = self.pbc - - D, D_len = self.get_distances(R, cell=cell, pbc=pbc) - - if vector: - return D - else: - return D_len - - def set_distance( - self, - a0, - a1, - distance, - fix=0.5, - mic=False, - mask=None, - indices=None, - add=False, - factor=False, - ): - """Set the distance between two atoms. - - Set the distance between atoms *a0* and *a1* to *distance*. - By default, the center of the two atoms will be fixed. Use - *fix=0* to fix the first atom, *fix=1* to fix the second - atom and *fix=0.5* (default) to fix the center of the bond. - - If *mask* or *indices* are set (*mask* overwrites *indices*), - only the atoms defined there are moved - (see :meth:`ase.Atoms.set_dihedral`). - - When *add* is true, the distance is changed by the value given. - In combination - with *factor* True, the value given is a factor scaling the distance. - - It is assumed that the atoms in *mask*/*indices* move together - with *a1*. If *fix=1*, only *a0* will therefore be moved.""" - - if a0 % len(self) == a1 % len(self): - raise ValueError('a0 and a1 must not be the same') - - if add: - oldDist = self.get_distance(a0, a1, mic=mic) - if factor: - newDist = oldDist * distance - else: - newDist = oldDist + distance - self.set_distance( - a0, - a1, - newDist, - fix=fix, - mic=mic, - mask=mask, - indices=indices, - add=False, - factor=False, - ) - return - - R = self.arrays['positions'] - D = np.array([R[a1] - R[a0]]) - - if mic: - D, D_len = self.find_mic(D, self.cell, self.pbc) + if wrap: + if 'pbc' not in wrap_kw: + wrap_kw['pbc'] = self.pbc + return ase.geometry.wrap_positions(self.positions, self.cell, **wrap_kw) else: - D_len = np.array([np.sqrt((D**2).sum())]) - x = 1.0 - distance / D_len[0] - - if mask is None and indices is None: - indices = [a0, a1] - elif mask: - indices = [i for i in range(len(self)) if mask[i]] - - for i in indices: - if i == a0: - R[a0] += (x * fix) * D[0] - else: - R[i] -= (x * (1.0 - fix)) * D[0] + return self.arrays['positions'].copy() def get_scaled_positions(self, wrap=True): """Get positions relative to unit cell. @@ -1633,17 +435,6 @@ def wrap(self, **wrap_kw): self.positions[:] = self.get_positions(wrap=True, **wrap_kw) - # @deprecated('Please use atoms.cell.volume') - # We kind of want to deprecate this, but the ValueError behaviour - # might be desirable. Should we do this? - def get_volume(self): - """Get volume of unit cell.""" - if self.cell.rank != 3: - raise ValueError( - f'You have {self.cell.rank} lattice vectors: volume not defined' - ) - return self.cell.volume - def _get_positions(self): """Return reference to positions-array for in-place manipulations.""" return self.arrays['positions'] From ed213717eeb27d3b3e4afffb4f621b82840b7b26 Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Wed, 16 Oct 2024 18:58:48 +0200 Subject: [PATCH 07/16] Converted lattice vectors to Nomad-compatible matrix format. --- .../schema_packages/model_system.py | 1 + .../schema_packages/particles_state.py | 33 ++----------------- 2 files changed, 4 insertions(+), 30 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 75205d95..259f751f 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -600,6 +600,7 @@ def to_particles(self, logger: 'BoundLogger') -> Optional[Particles]: # Lattice vectors if self.lattice_vectors is not None: particles.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + print(particles.get_cell()) else: logger.info('Could not find `ParticleCell.lattice_vectors`.') diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index b4a9658a..28afb1e6 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -282,47 +282,20 @@ def set_cell(self, cell, scale_atoms=False, apply_constraint=True): Parameters: - cell: 3x3 matrix or length 3 or 6 vector - Unit cell. A 3x3 matrix (the three unit cell vectors) or - just three numbers for an orthorhombic cell. Another option is - 6 numbers, which describes unit cell with lengths of unit cell - vectors and with angles between them (in degrees), in following - order: [len(a), len(b), len(c), angle(b,c), angle(a,c), - angle(a,b)]. First vector will lie in x-direction, second in + cell: 3x3 matrix + Unit cell. A 3x3 matrix (the three unit cell vectors). + First vector will lie in x-direction, second in xy-plane, and the third one in z-positive subspace. scale_atoms: bool Fix atomic positions or move atoms with the unit cell? Default behavior is to *not* move the atoms (scale_atoms=False). apply_constraint: bool Whether to apply constraints to the given cell. - - Examples: - - Two equivalent ways to define an orthorhombic cell: - - >>> atoms = Atoms('He') - >>> a, b, c = 7, 7.5, 8 - >>> atoms.set_cell([a, b, c]) - >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)]) - - FCC unit cell: - - >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)]) - - Hexagonal unit cell: - - >>> atoms.set_cell([a, a, c, 90, 90, 120]) - - Rhombohedral unit cell: - - >>> alpha = 77 - >>> atoms.set_cell([a, a, a, alpha, alpha, alpha]) """ # Override pbcs if and only if given a Cell object: cell = ase.Cell.new(cell) - # XXX not working well during initialize due to missing _constraints if apply_constraint and hasattr(self, '_constraints'): for constraint in self.constraints: if hasattr(constraint, 'adjust_cell'): From da38c467f901e6f8ba8702bb5df4db4df5e3f8bd Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Thu, 24 Oct 2024 18:26:37 +0200 Subject: [PATCH 08/16] Fixed wrong shape of bond_list in model_system --- src/nomad_simulations/schema_packages/model_system.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 259f751f..3725ebe4 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -1180,6 +1180,7 @@ class ModelSystem(System): # TODO improve description and add an example using the case in atom_indices bond_list = Quantity( type=np.int32, + shape=['*', 2], description=""" List of pairs of atom indices corresponding to bonds (e.g., as defined by a force field) within this atoms_group. From 4f1692ab806ab5f9762db68cbe14b6d11af4f9da Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Fri, 22 Nov 2024 13:57:37 +0100 Subject: [PATCH 09/16] Cell is 'ParticleCell', normalizer still calls atomic_cell.to_ase_atoms --- .../schema_packages/model_system.py | 5 +++- .../schema_packages/particles_state.py | 24 ------------------- 2 files changed, 4 insertions(+), 25 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 3725ebe4..68173b6f 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -521,7 +521,7 @@ class ParticleCell(Cell): def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs): super().__init__(m_def, m_context, **kwargs) # Set the name of the section - self.name = self.m_def.name + self.name = self.m_def.name #! self.name here is 'ParticleCell' def is_equal_cell(self, other) -> bool: """ @@ -1146,6 +1146,7 @@ class ModelSystem(System): ) cell = SubSection(sub_section=Cell.m_def, repeats=True) + print(f'cell: {Cell.m_def}') symmetry = SubSection(sub_section=Symmetry.m_def, repeats=True) @@ -1279,6 +1280,8 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: 'Could not find the originally parsed atomic system. `Symmetry` and `ChemicalFormula` extraction is thus not run.' ) return + + #! self.cell[0].name is ParticleCell if self.cell[0].name == 'AtomicCell': self.cell[0].type = 'original' ase_atoms = self.cell[0].to_ase_atoms(logger=logger) diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index 28afb1e6..acc419bb 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -94,30 +94,6 @@ class Particles: Periodic boundary conditions flags. Examples: True, False, 0, 1, (1, 1, 0), (True, False, False). Default value: False. - - Examples: - - These three are equivalent: - - >>> d = 1.104 # N2 bondlength - >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)]) - >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)]) - >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))]) - - FCC gold: - - >>> a = 4.05 # Gold lattice constant - >>> b = a / 2 - >>> fcc = Atoms('Au', - ... cell=[(0, b, b), (b, 0, b), (b, b, 0)], - ... pbc=True) - - Hydrogen wire: - - >>> d = 0.9 # H-H distance - >>> h = Atoms('H', positions=[(0, 0, 0)], - ... cell=(d, 0, 0), - ... pbc=(1, 0, 0)) """ def __init__( From a3805c214ca0dab7a7c3dd5739b2dd457863801a Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Wed, 27 Nov 2024 16:08:21 +0100 Subject: [PATCH 10/16] Make chemical formula generation and normalization of atom types conditional --- .../schema_packages/model_system.py | 20 ++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 68173b6f..5363b6df 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -1,4 +1,5 @@ import re +import sys from typing import TYPE_CHECKING, Optional import ase @@ -1304,11 +1305,16 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: sec_symmetry = self.m_create(Symmetry) sec_symmetry.normalize(archive, logger) + print(self.model_system) + print(self.cell[0]) + print(self.cell[0].particles_state[0].particle_type) + sys.exit() # Creating and normalizing ChemicalFormula section - # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) - sec_chemical_formula = self.m_create(ChemicalFormula) - sec_chemical_formula.normalize(archive, logger) - if sec_chemical_formula.m_cache: - self.elemental_composition = sec_chemical_formula.m_cache.get( - 'elemental_composition', [] - ) + if self.cell[0].name == 'AtomicCell': + # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) + sec_chemical_formula = self.m_create(ChemicalFormula) + sec_chemical_formula.normalize(archive, logger) + if sec_chemical_formula.m_cache: + self.elemental_composition = sec_chemical_formula.m_cache.get( + 'elemental_composition', [] + ) From 2bd3449352520bdfe5067582a16dadfea57aa110 Mon Sep 17 00:00:00 2001 From: ndaelman-hu <107392603+ndaelman-hu@users.noreply.github.com> Date: Wed, 16 Oct 2024 14:36:00 +0200 Subject: [PATCH 11/16] Resolve conflicts during rebase to "cg-particle-support" --------- Co-authored-by: ndaelman --- .../schema_packages/__init__.py | 4 +- .../schema_packages/model_system.py | 208 +++++++++++------ .../schema_packages/utils/__init__.py | 1 + .../schema_packages/utils/utils.py | 18 +- tests/test_model_system.py | 210 +++++++++++------- 5 files changed, 281 insertions(+), 160 deletions(-) diff --git a/src/nomad_simulations/schema_packages/__init__.py b/src/nomad_simulations/schema_packages/__init__.py index 8b730793..78d66557 100644 --- a/src/nomad_simulations/schema_packages/__init__.py +++ b/src/nomad_simulations/schema_packages/__init__.py @@ -31,8 +31,8 @@ class NOMADSimulationsEntryPoint(SchemaPackageEntryPoint): description='Limite of the number of atoms in the unit cell to be treated for the system type classification from MatID to work. This is done to avoid overhead of the package.', ) equal_cell_positions_tolerance: float = Field( - 1e-12, - description='Tolerance (in meters) for the cell positions to be considered equal.', + 12, + description='Decimal order or tolerance (in meters) for comparing cell positions.', ) def load(self): diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 5363b6df..0c307d42 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -1,5 +1,24 @@ +# +# Copyright The NOMAD Authors. +# +# This file is part of NOMAD. See https://nomad-lab.eu for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + import re -import sys +from functools import lru_cache +from hashlib import sha1 from typing import TYPE_CHECKING, Optional import ase @@ -23,6 +42,10 @@ from nomad.units import ureg if TYPE_CHECKING: + from collections.abc import Generator + from typing import Any, Callable, Optional + + import pint from nomad.datamodel.datamodel import EntryArchive from nomad.metainfo import Context, Section from structlog.stdlib import BoundLogger @@ -30,6 +53,7 @@ from nomad_simulations.schema_packages.atoms_state import AtomsState from nomad_simulations.schema_packages.particles_state import Particles, ParticlesState from nomad_simulations.schema_packages.utils import ( + catch_not_implemented, get_sibling_section, is_not_representative, ) @@ -202,6 +226,72 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: return +def _check_implemented(func: 'Callable'): + """ + Decorator to restrict the comparison functions to the same class. + """ + + def wrapper(self, other): + if not isinstance(other, self.__class__): + return NotImplemented + return func(self, other) + + return wrapper + + +class PartialOrderElement: + def __init__(self, representative_variable): + self.representative_variable = representative_variable + + def __hash__(self): + return self.representative_variable.__hash__() + + @_check_implemented + def __eq__(self, other): + return self.representative_variable == other.representative_variable + + @_check_implemented + def __lt__(self, other): + return False + + @_check_implemented + def __gt__(self, other): + return False + + def __le__(self, other): + return self.__eq__(other) + + def __ge__(self, other): + return self.__eq__(other) + + # __ne__ assumes that usage in a finite set with its comparison definitions + + +class HashedPositions(PartialOrderElement): + # `representative_variable` is a `pint.Quantity` object + + def __hash__(self): + hash_str = sha1( + np.ascontiguousarray( + np.round( + self.representative_variable.to_base_units().magnitude, + decimals=configuration.equal_cell_positions_tolerance, + out=None, + ) + ).tobytes() + ).hexdigest() + return int(hash_str, 16) + + def __eq__(self, other): + """Equality as defined between HashedPositions.""" + if ( + self.representative_variable is None + or other.representative_variable is None + ): + return NotImplemented + return np.allclose(self.representative_variable, other.representative_variable) + + class Cell(GeometricSpace): """ A base section used to specify the cell quantities of a system at a given moment in time. @@ -219,7 +309,7 @@ class Cell(GeometricSpace): type=MEnum('original', 'primitive', 'conventional'), description=""" Representation type of the cell structure. It might be: - - 'original' as in origanally parsed, + - 'original' as in originally parsed, - 'primitive' as the primitive unit cell, - 'conventional' as the conventional cell used for referencing. """, @@ -280,45 +370,36 @@ class Cell(GeometricSpace): """, ) - def _check_positions(self, positions_1, positions_2) -> list: - # Check that all the `positions`` of `cell_1` match with the ones in `cell_2` - check_positions = [] - for i1, pos1 in enumerate(positions_1): - for i2, pos2 in enumerate(positions_2): - if np.allclose( - pos1, pos2, atol=configuration.equal_cell_positions_tolerance - ): - check_positions.append([i1, i2]) - break - return check_positions + @staticmethod + def _generate_comparer(obj: 'Cell') -> 'Generator[Any, None, None]': + try: + return ((HashedPositions(pos)) for pos in obj.positions) + except AttributeError: + raise NotImplementedError - def is_equal_cell(self, other) -> bool: - """ - Check if the cell is equal to an`other` cell by comparing the `positions`. - Args: - other: The other cell to compare with. - Returns: - bool: True if the cells are equal, False otherwise. - """ - # TODO implement checks on `lattice_vectors` and other quantities to ensure the equality of primitive cells - if not isinstance(other, Cell): - return False + @catch_not_implemented + def is_lt_cell(self, other) -> bool: + return set(self._generate_comparer(self)) < set(self._generate_comparer(other)) - # If the `positions` are empty, return False - if self.positions is None or other.positions is None: - return False + @catch_not_implemented + def is_gt_cell(self, other) -> bool: + return set(self._generate_comparer(self)) > set(self._generate_comparer(other)) - # The `positions` should have the same length (same number of positions) - if len(self.positions) != len(other.positions): - return False - n_positions = len(self.positions) + @catch_not_implemented + def is_le_cell(self, other) -> bool: + return set(self._generate_comparer(self)) <= set(self._generate_comparer(other)) - check_positions = self._check_positions( - positions_1=self.positions, positions_2=other.positions - ) - if len(check_positions) != n_positions: - return False - return True + @catch_not_implemented + def is_ge_cell(self, other) -> bool: + return set(self._generate_comparer(self)) >= set(self._generate_comparer(other)) + + @catch_not_implemented + def is_equal_cell(self, other) -> bool: # TODO: improve naming + return set(self._generate_comparer(self)) == set(self._generate_comparer(other)) + + def is_ne_cell(self, other) -> bool: + # this does not hold in general, but here we use finite sets + return not self.is_equal_cell(other) def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: super().normalize(archive, logger) @@ -363,40 +444,20 @@ def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwarg # Set the name of the section self.name = self.m_def.name - def is_equal_cell(self, other) -> bool: - """ - Check if the atomic cell is equal to an`other` atomic cell by comparing the `positions` and - the `AtomsState[*].chemical_symbol`. - Args: - other: The other atomic cell to compare with. - Returns: - bool: True if the atomic cells are equal, False otherwise. - """ - if not isinstance(other, AtomicCell): - return False - - # Compare positions using the parent sections's `__eq__` method - if not super().is_equal_cell(other=other): - return False - - # Check that the `chemical_symbol` of the atoms in `cell_1` match with the ones in `cell_2` - check_positions = self._check_positions( - positions_1=self.positions, positions_2=other.positions - ) + @staticmethod + def _generate_comparer(obj: 'AtomicCell') -> 'Generator[Any, None, None]': + # presumes `atoms_state` mapping 1-to-1 with `positions` and conserves the order try: - for atom in check_positions: - element_1 = self.atoms_state[atom[0]].chemical_symbol - element_2 = other.atoms_state[atom[1]].chemical_symbol - if element_1 != element_2: - return False - except Exception: - return False - return True + return ( + (HashedPositions(pos), PartialOrderElement(st.chemical_symbol)) + for pos, st in zip(obj.positions, obj.atoms_state) + ) + except AttributeError: + raise NotImplementedError def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]: """ Get the chemical symbols of the atoms in the atomic cell. These are defined on `atoms_state[*].chemical_symbol`. - Args: logger (BoundLogger): The logger to log messages. @@ -414,7 +475,7 @@ def get_chemical_symbols(self, logger: 'BoundLogger') -> list[str]: chemical_symbols.append(atom_state.chemical_symbol) return chemical_symbols - def to_ase_atoms(self, logger: 'BoundLogger') -> Optional[ase.Atoms]: + def to_ase_atoms(self, logger: 'BoundLogger') -> 'Optional[ase.Atoms]': """ Generates an ASE Atoms object with the most basic information from the parsed `AtomicCell` section (labels, periodic_boundary_conditions, positions, and lattice_vectors). @@ -766,8 +827,11 @@ class Symmetry(ArchiveSection): ) def resolve_analyzed_atomic_cell( - self, symmetry_analyzer: SymmetryAnalyzer, cell_type: str, logger: 'BoundLogger' - ) -> Optional[AtomicCell]: + self, + symmetry_analyzer: 'SymmetryAnalyzer', + cell_type: str, + logger: 'BoundLogger', + ) -> 'Optional[AtomicCell]': """ Resolves the `AtomicCell` section from the `SymmetryAnalyzer` object and the cell_type (primitive or conventional). @@ -811,8 +875,8 @@ def resolve_analyzed_atomic_cell( return atomic_cell def resolve_bulk_symmetry( - self, original_atomic_cell: AtomicCell, logger: 'BoundLogger' - ) -> tuple[Optional[AtomicCell], Optional[AtomicCell]]: + self, original_atomic_cell: 'AtomicCell', logger: 'BoundLogger' + ) -> 'tuple[Optional[AtomicCell], Optional[AtomicCell]]': """ Resolves the symmetry of the material being simulated using MatID and the originally parsed data under original_atomic_cell. It generates two other @@ -1308,7 +1372,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: print(self.model_system) print(self.cell[0]) print(self.cell[0].particles_state[0].particle_type) - sys.exit() + # sys.exit() # Creating and normalizing ChemicalFormula section if self.cell[0].name == 'AtomicCell': # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) diff --git a/src/nomad_simulations/schema_packages/utils/__init__.py b/src/nomad_simulations/schema_packages/utils/__init__.py index 52d9ca22..f9945a34 100644 --- a/src/nomad_simulations/schema_packages/utils/__init__.py +++ b/src/nomad_simulations/schema_packages/utils/__init__.py @@ -1,5 +1,6 @@ from .utils import ( RussellSaundersState, + catch_not_implemented, get_composition, get_sibling_section, get_variables, diff --git a/src/nomad_simulations/schema_packages/utils/utils.py b/src/nomad_simulations/schema_packages/utils/utils.py index 1d40aa4a..eff18376 100644 --- a/src/nomad_simulations/schema_packages/utils/utils.py +++ b/src/nomad_simulations/schema_packages/utils/utils.py @@ -5,7 +5,7 @@ from nomad.config import config if TYPE_CHECKING: - from typing import Optional + from typing import Callable, Optional from nomad.datamodel.data import ArchiveSection from structlog.stdlib import BoundLogger @@ -154,3 +154,19 @@ def get_composition(children_names: 'list[str]') -> str: children_count_tup = np.unique(children_names, return_counts=True) formula = ''.join([f'{name}({count})' for name, count in zip(*children_count_tup)]) return formula if formula else None + + +def catch_not_implemented(func: 'Callable') -> 'Callable': + """ + Decorator to default comparison functions outside the same class to `False`. + """ + + def wrapper(self, other) -> bool: + if not isinstance(other, self.__class__): + return False # ? should this throw an error instead? + try: + return func(self, other) + except (TypeError, NotImplementedError): + return False + + return wrapper diff --git a/tests/test_model_system.py b/tests/test_model_system.py index f334da23..088ecc6b 100644 --- a/tests/test_model_system.py +++ b/tests/test_model_system.py @@ -18,96 +18,104 @@ from .conftest import generate_atomic_cell -class TestCell: +class TestAtomicCell: """ - Test the `Cell` section defined in model_system.py + Test the `AtomicCell`, `Cell` and `GeometricSpace` classes defined in model_system.py """ @pytest.mark.parametrize( 'cell_1, cell_2, result', [ - (Cell(), None, False), # one cell is None - (Cell(), Cell(), False), # both cells are empty + (Cell(), None, {'lt': False, 'gt': False, 'eq': False}), # one cell is None + # (Cell(), Cell(), False), # both cells are empty + # ( + # Cell(positions=[[1, 0, 0]]), + # Cell(), + # False, + # ), # one cell has positions, the other is empty ( Cell(positions=[[1, 0, 0]]), - Cell(), - False, - ), # one cell has positions, the other is empty + Cell(positions=[[2, 0, 0]]), + {'lt': False, 'gt': False, 'eq': False}, + ), # position vectors are treated as the fundamental set elements ( Cell(positions=[[1, 0, 0], [0, 1, 0]]), Cell(positions=[[1, 0, 0]]), - False, - ), # length mismatch - ( - Cell(positions=[[1, 0, 0], [0, 1, 0]]), - Cell(positions=[[1, 0, 0], [0, -1, 0]]), - False, - ), # different positions - ( - Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - True, - ), # same ordered positions - ( - Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - Cell(positions=[[1, 0, 0], [0, 0, 1], [0, 1, 0]]), - True, - ), # different ordered positions but same cell - ], - ) - def test_is_equal_cell(self, cell_1: Cell, cell_2: Cell, result: bool): - """ - Test the `is_equal_cell` methods of `Cell`. - """ - assert cell_1.is_equal_cell(other=cell_2) == result - - -class TestAtomicCell: - """ - Test the `AtomicCell`, `Cell` and `GeometricSpace` classes defined in model_system.py - """ - - @pytest.mark.parametrize( - 'cell_1, cell_2, result', - [ - (Cell(), None, False), # one cell is None - (Cell(), Cell(), False), # both cells are empty + {'lt': False, 'gt': True, 'eq': False}, + ), # one is a subset of the other ( Cell(positions=[[1, 0, 0]]), - Cell(), - False, - ), # one cell has positions, the other is empty - ( Cell(positions=[[1, 0, 0], [0, 1, 0]]), - Cell(positions=[[1, 0, 0]]), - False, - ), # length mismatch + {'lt': True, 'gt': False, 'eq': False}, + ), # one is a subset of the other ( Cell(positions=[[1, 0, 0], [0, 1, 0]]), Cell(positions=[[1, 0, 0], [0, -1, 0]]), - False, + {'lt': False, 'gt': False, 'eq': False}, ), # different positions ( Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - True, + {'lt': False, 'gt': False, 'eq': True}, ), # same ordered positions ( Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), Cell(positions=[[1, 0, 0], [0, 0, 1], [0, 1, 0]]), - True, + {'lt': False, 'gt': False, 'eq': True}, ), # different ordered positions but same cell + # ( + # AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), + # Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), + # False, + # ), # one atomic cell and another cell (missing chemical symbols) + # ( + # AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), + # AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), + # False, + # ), # missing chemical symbols + # ND: the comparison will now return an error here + # handling a case that should be resolved by the normalizer falls outside its scope ( - AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - Cell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - False, - ), # one atomic cell and another cell (missing chemical symbols) + AtomicCell( + positions=[[1, 0, 0]], + atoms_state=[ + AtomsState(chemical_symbol='O'), + ], + ), + AtomicCell( + positions=[[1, 0, 0]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + ], + ), + {'lt': False, 'gt': False, 'eq': False}, + ), # chemical symbols are treated as the fundamental set elements ( - AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - AtomicCell(positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]]), - False, - ), # missing chemical symbols + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='O'), + ], + ), + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + ], + ), + {'lt': False, 'gt': True, 'eq': False}, + ), # one is a subset of the other ( + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + ], + ), AtomicCell( positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], atoms_state=[ @@ -116,6 +124,16 @@ class TestAtomicCell: AtomsState(chemical_symbol='O'), ], ), + {'lt': True, 'gt': False, 'eq': False}, + ), # one is a subset of the other + ( + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='O'), + ], + ), AtomicCell( positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], atoms_state=[ @@ -124,7 +142,26 @@ class TestAtomicCell: AtomsState(chemical_symbol='O'), ], ), - True, + {'lt': False, 'gt': False, 'eq': False}, + ), + ( + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='O'), + ], + ), + AtomicCell( + positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='O'), + ], + ), + {'lt': False, 'gt': False, 'eq': True}, ), # same ordered positions and chemical symbols ( AtomicCell( @@ -143,7 +180,7 @@ class TestAtomicCell: AtomsState(chemical_symbol='O'), ], ), - False, + {'lt': False, 'gt': False, 'eq': False}, ), # same ordered positions but different chemical symbols ( AtomicCell( @@ -162,38 +199,41 @@ class TestAtomicCell: AtomsState(chemical_symbol='H'), ], ), - True, - ), # different ordered positions but same chemical symbols - ], - ) - def test_is_equal_cell(self, cell_1: Cell, cell_2: Cell, result: bool): - """ - Test the `is_equal_cell` methods of `AtomicCell`. - """ - assert cell_1.is_equal_cell(other=cell_2) == result - - @pytest.mark.parametrize( - 'atomic_cell, result', - [ - (AtomicCell(), []), - (AtomicCell(atoms_state=[AtomsState(chemical_symbol='H')]), ['H']), + {'lt': False, 'gt': False, 'eq': True}, + ), # same position-symbol map, different overall order ( AtomicCell( + positions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], + atoms_state=[ + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='H'), + AtomsState(chemical_symbol='O'), + ], + ), + AtomicCell( + positions=[[1, 0, 0], [0, 0, 1], [0, 1, 0]], atoms_state=[ AtomsState(chemical_symbol='H'), - AtomsState(chemical_symbol='Fe'), + AtomsState(chemical_symbol='H'), AtomsState(chemical_symbol='O'), - ] + ], ), - ['H', 'Fe', 'O'], - ), + {'lt': False, 'gt': False, 'eq': False}, + ), # different position-symbol map ], ) - def test_get_chemical_symbols(self, atomic_cell: AtomicCell, result: list[str]): + def test_partial_order( + self, cell_1: 'Cell', cell_2: 'Cell', result: dict[str, bool] + ): """ - Test the `get_chemical_symbols` method of `AtomicCell`. + Test the comparison operators of `Cell` and `AtomicCell`. """ - assert atomic_cell.get_chemical_symbols(logger=logger) == result + assert cell_1.is_lt_cell(cell_2) == result['lt'] + assert cell_1.is_gt_cell(cell_2) == result['gt'] + assert cell_1.is_le_cell(cell_2) == (result['lt'] or result['eq']) + assert cell_1.is_ge_cell(cell_2) == (result['gt'] or result['eq']) + assert cell_1.is_equal_cell(cell_2) == result['eq'] + assert cell_1.is_ne_cell(cell_2) == (not result['eq']) @pytest.mark.parametrize( 'chemical_symbols, atomic_numbers, formula, lattice_vectors, positions, periodic_boundary_conditions', From 763e5e9631f74607c05f655a3553a63f2af29721 Mon Sep 17 00:00:00 2001 From: ndaelman-hu <107392603+ndaelman-hu@users.noreply.github.com> Date: Mon, 18 Nov 2024 11:06:50 +0100 Subject: [PATCH 12/16] Bump dependencies (#147) Co-authored-by: ndaelman --- pyproject.toml | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0cd48773..ae813b12 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,12 +45,13 @@ dependencies = [ [project.optional-dependencies] dev = [ - "mypy==1.0.1", - "ruff", - "pytest", - "pytest-timeout", - "pytest-cov", - "structlog", + 'mypy==1.0.1', + 'pytest>= 5.3.0, <8', + 'pytest-timeout>=1.4.2', + 'pytest-cov>=2.7.1', + 'ruff>=0.6', + 'structlog>=1.0', + 'typing-extensions>=4.12', ] [tool.uv] From 91fe83b23c20ad55e8b8a82c6f3d4029b77c359e Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Fri, 11 Oct 2024 15:28:21 +0200 Subject: [PATCH 13/16] fixing rebase --- .../schema_packages/model_system.py | 87 +++-- .../schema_packages/particles_state.py | 329 ++++++++++++++---- 2 files changed, 296 insertions(+), 120 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 0c307d42..ce065f36 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -557,7 +557,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: class ParticleCell(Cell): """ - A base section used to specify the atomic cell information of a system. + A base section used to specify the particle cell information of a system. """ particles_state = SubSection(sub_section=ParticlesState.m_def, repeats=True) @@ -565,65 +565,66 @@ class ParticleCell(Cell): n_particles = Quantity( type=np.int32, description=""" - Number of atoms in the atomic cell. + Number of particles in the particle cell. """, ) equivalent_particles = Quantity( type=np.int32, - shape=['n_atoms'], + shape=['n_particle'], description=""" - List of equivalent atoms as defined in `atoms`. If no equivalent atoms are found, - then the list is simply the index of each element, e.g.: - - [0, 1, 2, 3] all four atoms are non-equivalent. - - [0, 0, 0, 3] three equivalent atoms and one non-equivalent. + List of equivalent particles as defined in `particles`. If no equivalent particles + are found, then the list is simply the index of each element, e.g.: + - [0, 1, 2, 3] all four particles are non-equivalent. + - [0, 0, 0, 3] three equivalent particles and one non-equivalent. """, ) def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs): super().__init__(m_def, m_context, **kwargs) # Set the name of the section - self.name = self.m_def.name #! self.name here is 'ParticleCell' - - def is_equal_cell(self, other) -> bool: - """ - Check if the atomic cell is equal to an`other` atomic cell by comparing the `positions` and - the `AtomsState[*].chemical_symbol`. - Args: - other: The other atomic cell to compare with. - Returns: - bool: True if the atomic cells are equal, False otherwise. - """ - if not isinstance(other, ParticleCell): - return False - - # Compare positions using the parent sections's `__eq__` method - if not super().is_equal_cell(other=other): - return False + self.name = self.m_def.name - # Check that the `chemical_symbol` of the atoms in `cell_1` match with the ones in `cell_2` - check_positions = self._check_positions( - positions_1=self.positions, positions_2=other.positions - ) - try: - for particle in check_positions: - type_1 = self.particles_state[particle[0]].particle_type - type_2 = other.particles_state[particle[1]].particle_type - if type_1 != type_2: - return False - except Exception: - return False - return True + # def is_equal_cell(self, other) -> bool: + # """ + # Check if the atomic cell is equal to an`other` atomic cell by comparing the `positions` and + # the `AtomsState[*].chemical_symbol`. + # Args: + # other: The other atomic cell to compare with. + # Returns: + # bool: True if the atomic cells are equal, False otherwise. + # """ + # if not isinstance(other, ParticleCell): + # return False + + # # Compare positions using the parent sections's `__eq__` method + # if not super().is_equal_cell(other=other): + # return False + + # # Check that the `chemical_symbol` of the atoms in `cell_1` match with the ones in `cell_2` + # check_positions = self._check_positions( + # positions_1=self.positions, positions_2=other.positions + # ) + # try: + # for particle in check_positions: + # type_1 = self.particles_state[particle[0]].particle_type + # type_2 = other.particles_state[particle[1]].particle_type + # if type_1 != type_2: + # return False + # except Exception: + # return False + # return True def get_particle_types(self, logger: 'BoundLogger') -> list[str]: """ - Get the chemical symbols of the atoms in the atomic cell. These are defined on `atoms_state[*].chemical_symbol`. + Get the chemical symbols of the particle in the particle cell. + These are defined on `particles_state[*].chemical_symbol`. Args: logger (BoundLogger): The logger to log messages. Returns: - list: The list of chemical symbols of the atoms in the atomic cell. + list: The list of chemical symbols of the particles in the particle cell. """ if not self.particles_state: return [] @@ -1346,7 +1347,6 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: ) return - #! self.cell[0].name is ParticleCell if self.cell[0].name == 'AtomicCell': self.cell[0].type = 'original' ase_atoms = self.cell[0].to_ase_atoms(logger=logger) @@ -1369,12 +1369,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: sec_symmetry = self.m_create(Symmetry) sec_symmetry.normalize(archive, logger) - print(self.model_system) - print(self.cell[0]) - print(self.cell[0].particles_state[0].particle_type) - # sys.exit() - # Creating and normalizing ChemicalFormula section - if self.cell[0].name == 'AtomicCell': + # Creating and normalizing ChemicalFormula section # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) sec_chemical_formula = self.m_create(ChemicalFormula) sec_chemical_formula.normalize(archive, logger) diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index acc419bb..f3504708 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -94,6 +94,30 @@ class Particles: Periodic boundary conditions flags. Examples: True, False, 0, 1, (1, 1, 0), (True, False, False). Default value: False. + + Examples: + + These three are equivalent: + + >>> d = 1.104 # N2 bondlength + >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)]) + >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)]) + >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))]) + + FCC gold: + + >>> a = 4.05 # Gold lattice constant + >>> b = a / 2 + >>> fcc = Atoms('Au', + ... cell=[(0, b, b), (b, 0, b), (b, b, 0)], + ... pbc=True) + + Hydrogen wire: + + >>> d = 0.9 # H-H distance + >>> h = Atoms('H', positions=[(0, 0, 0)], + ... cell=(d, 0, 0), + ... pbc=(1, 0, 0)) """ def __init__( @@ -120,79 +144,79 @@ def __init__( particles = None - # if hasattr(types, 'get_positions'): - # atoms = types - # types = None - # elif ( - # isinstance(types, (list, tuple)) - # and len(types) > 0 - # and isinstance(types[0], Atom) - # ): - # # Get data from a list or tuple of Atom objects: - # data = [ - # [atom.get_raw(name) for atom in types] - # for name in [ - # 'position', - # 'number', - # 'tag', - # 'momentum', - # 'mass', - # 'magmom', - # 'charge', - # ] - # ] - # atoms = self.__class__(None, *data) - # types = None - - # if atoms is not None: - # # Get data from another Atoms object: - # if scaled_positions is not None: - # raise NotImplementedError - # if types is None and numbers is None: - # numbers = atoms.get_atomic_numbers() - # if positions is None: - # positions = atoms.get_positions() - # if tags is None and atoms.has('tags'): - # tags = atoms.get_tags() - # if momenta is None and atoms.has('momenta'): - # momenta = atoms.get_momenta() - # if magmoms is None and atoms.has('initial_magmoms'): - # magmoms = atoms.get_initial_magnetic_moments() - # if masses is None and atoms.has('masses'): - # masses = atoms.get_masses() - # if charges is None and atoms.has('initial_charges'): - # charges = atoms.get_initial_charges() - # if cell is None: - # cell = atoms.get_cell() - # if celldisp is None: - # celldisp = atoms.get_celldisp() - # if pbc is None: - # pbc = atoms.get_pbc() - - # self.arrays = {} - - # if types is None: - # if numbers is None: - # if positions is not None: - # natoms = len(positions) - # elif scaled_positions is not None: - # natoms = len(scaled_positions) - # else: - # natoms = 0 - # numbers = np.zeros(natoms, int) - # self.new_array('numbers', numbers, int) - # else: - # if numbers is not None: - # raise TypeError('Use only one of "types" and "numbers".') - # else: - # self.new_array('numbers', types2numbers(types), int) - - # if self.numbers.ndim != 1: - # raise ValueError('"numbers" must be 1-dimensional.') - - if cell is None: - cell = np.zeros((3, 3)) - self.set_cell(cell) + # if hasattr(symbols, 'get_positions'): + # atoms = symbols + # symbols = None + # elif ( + # isinstance(symbols, (list, tuple)) + # and len(symbols) > 0 + # and isinstance(symbols[0], Atom) + # ): + # # Get data from a list or tuple of Atom objects: + # data = [ + # [atom.get_raw(name) for atom in symbols] + # for name in [ + # 'position', + # 'number', + # 'tag', + # 'momentum', + # 'mass', + # 'magmom', + # 'charge', + # ] + # ] + # atoms = self.__class__(None, *data) + # symbols = None + + # if atoms is not None: + # # Get data from another Atoms object: + # if scaled_positions is not None: + # raise NotImplementedError + # if symbols is None and numbers is None: + # numbers = atoms.get_atomic_numbers() + # if positions is None: + # positions = atoms.get_positions() + # if tags is None and atoms.has('tags'): + # tags = atoms.get_tags() + # if momenta is None and atoms.has('momenta'): + # momenta = atoms.get_momenta() + # if magmoms is None and atoms.has('initial_magmoms'): + # magmoms = atoms.get_initial_magnetic_moments() + # if masses is None and atoms.has('masses'): + # masses = atoms.get_masses() + # if charges is None and atoms.has('initial_charges'): + # charges = atoms.get_initial_charges() + # if cell is None: + # cell = atoms.get_cell() + # if celldisp is None: + # celldisp = atoms.get_celldisp() + # if pbc is None: + # pbc = atoms.get_pbc() + + # self.arrays = {} + + # if symbols is None: + # if numbers is None: + # if positions is not None: + # natoms = len(positions) + # elif scaled_positions is not None: + # natoms = len(scaled_positions) + # else: + # natoms = 0 + # numbers = np.zeros(natoms, int) + # self.new_array('numbers', numbers, int) + # else: + # if numbers is not None: + # raise TypeError('Use only one of "symbols" and "numbers".') + # else: + # self.new_array('numbers', symbols2numbers(symbols), int) + + # if self.numbers.ndim != 1: + # raise ValueError('"numbers" must be 1-dimensional.') + + # if cell is None: + # cell = np.zeros((3, 3)) + # self.set_cell(cell) # if celldisp is None: # celldisp = np.zeros(shape=(3, 1)) @@ -206,7 +230,7 @@ def __init__( # positions = np.dot(scaled_positions, self.cell) # else: # if scaled_positions is not None: - # raise TypeError('Use only one of "types" and "numbers".') + # raise TypeError('Use only one of "symbols" and "numbers".') # self.new_array('positions', positions, float, (3,)) # self.set_tags(default(tags, 0)) # self.set_masses(default(masses, None)) @@ -312,6 +336,152 @@ def get_pbc(self): """Get periodic boundary condition flags.""" return self.pbc.copy() + def new_array(self, name, a, dtype=None, shape=None): + """Add new array. + + If *shape* is not *None*, the shape of *a* will be checked.""" + + if dtype is not None: + a = np.array(a, dtype, order='C') + if len(a) == 0 and shape is not None: + a.shape = (-1,) + shape + else: + if not a.flags['C_CONTIGUOUS']: + a = np.ascontiguousarray(a) + else: + a = a.copy() + + if name in self.arrays: + raise RuntimeError(f'Array {name} already present') + + for b in self.arrays.values(): + if len(a) != len(b): + raise ValueError( + 'Array "%s" has wrong length: %d != %d.' % (name, len(a), len(b)) + ) + break + + if shape is not None and a.shape[1:] != shape: + raise ValueError( + 'Array "%s" has wrong shape %s != %s.' + % (name, a.shape, (a.shape[0:1] + shape)) + ) + + self.arrays[name] = a + + def get_array(self, name, copy=True): + """Get an array. + + Returns a copy unless the optional argument copy is false. + """ + if copy: + return self.arrays[name].copy() + else: + return self.arrays[name] + + def set_array(self, name, a, dtype=None, shape=None): + """Update array. + + If *shape* is not *None*, the shape of *a* will be checked. + If *a* is *None*, then the array is deleted.""" + + b = self.arrays.get(name) + if b is None: + if a is not None: + self.new_array(name, a, dtype, shape) + else: + if a is None: + del self.arrays[name] + else: + a = np.asarray(a) + if a.shape != b.shape: + raise ValueError( + 'Array "%s" has wrong shape %s != %s.' + % (name, a.shape, b.shape) + ) + b[:] = a + + def has(self, name): + """Check for existence of array. + + name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms', + 'initial_charges'.""" + # XXX extend has to calculator properties + return name in self.arrays + + def set_tags(self, tags): + """Set tags for all atoms. If only one tag is supplied, it is + applied to all atoms.""" + if isinstance(tags, int): + tags = [tags] * len(self) + self.set_array('tags', tags, int, ()) + + def get_tags(self): + """Get integer array of tags.""" + if 'tags' in self.arrays: + return self.arrays['tags'].copy() + else: + return np.zeros(len(self), int) + + def set_momenta(self, momenta, apply_constraint=True): + """Set momenta.""" + if apply_constraint and len(self.constraints) > 0 and momenta is not None: + momenta = np.array(momenta) # modify a copy + for constraint in self.constraints: + if hasattr(constraint, 'adjust_momenta'): + constraint.adjust_momenta(self, momenta) + self.set_array('momenta', momenta, float, (3,)) + + def set_velocities(self, velocities): + """Set the momenta by specifying the velocities.""" + self.set_momenta(self.get_masses()[:, np.newaxis] * velocities) + + def get_momenta(self): + """Get array of momenta.""" + if 'momenta' in self.arrays: + return self.arrays['momenta'].copy() + else: + return np.zeros((len(self), 3)) + + # def set_masses(self, masses='defaults'): + # """Set atomic masses in atomic mass units. + + # The array masses should contain a list of masses. In case + # the masses argument is not given or for those elements of the + # masses list that are None, standard values are set.""" + + # if isinstance(masses, str): + # if masses == 'defaults': + # masses = atomic_masses[self.arrays['numbers']] + # elif masses == 'most_common': + # masses = atomic_masses_common[self.arrays['numbers']] + # elif masses is None: + # pass + # elif not isinstance(masses, np.ndarray): + # masses = list(masses) + # for i, mass in enumerate(masses): + # if mass is None: + # masses[i] = atomic_masses[self.numbers[i]] + # self.set_array('masses', masses, float, ()) + + # def get_masses(self): + # """Get array of masses in atomic mass units.""" + # if 'masses' in self.arrays: + # return self.arrays['masses'].copy() + # else: + # return atomic_masses[self.arrays['numbers']] + + def get_charges(self): + """Get calculated charges.""" + if self._calc is None: + raise RuntimeError('Atoms object has no calculator.') + try: + return self._calc.get_charges(self) + except AttributeError: + from ase.calculators.calculator import PropertyNotImplementedError + + raise PropertyNotImplementedError + def set_positions(self, newpositions, apply_constraint=True): """Set positions, honoring any constraints. To ignore constraints, use *apply_constraint=False*.""" @@ -384,6 +554,17 @@ def wrap(self, **wrap_kw): self.positions[:] = self.get_positions(wrap=True, **wrap_kw) + # @deprecated('Please use atoms.cell.volume') + # We kind of want to deprecate this, but the ValueError behaviour + # might be desirable. Should we do this? + def get_volume(self): + """Get volume of unit cell.""" + if self.cell.rank != 3: + raise ValueError( + f'You have {self.cell.rank} lattice vectors: volume not defined' + ) + return self.cell.volume + def _get_positions(self): """Return reference to positions-array for in-place manipulations.""" return self.arrays['positions'] From 3604a86b4531cda1312f7280c12d3cacc8536b4b Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Mon, 14 Oct 2024 17:47:29 +0200 Subject: [PATCH 14/16] fixing rebase --- .../schema_packages/model_system.py | 13 +- .../schema_packages/particles_state.py | 305 +++++------------- 2 files changed, 80 insertions(+), 238 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index ce065f36..7f47f0bf 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -51,7 +51,7 @@ from structlog.stdlib import BoundLogger from nomad_simulations.schema_packages.atoms_state import AtomsState -from nomad_simulations.schema_packages.particles_state import Particles, ParticlesState +from nomad_simulations.schema_packages.particles_state import Particles, Particles, ParticlesState from nomad_simulations.schema_packages.utils import ( catch_not_implemented, get_sibling_section, @@ -660,12 +660,11 @@ def to_particles(self, logger: 'BoundLogger') -> Optional[Particles]: self.periodic_boundary_conditions = [False, False, False] particles.set_pbc(pbc=self.periodic_boundary_conditions) - # Lattice vectors - if self.lattice_vectors is not None: - particles.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) - print(particles.get_cell()) - else: - logger.info('Could not find `ParticleCell.lattice_vectors`.') + # # Lattice vectors + # if self.lattice_vectors is not None: + # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + # else: + # logger.info('Could not find `AtomicCell.lattice_vectors`.') # Positions if self.positions is not None: diff --git a/src/nomad_simulations/schema_packages/particles_state.py b/src/nomad_simulations/schema_packages/particles_state.py index f3504708..28afb1e6 100644 --- a/src/nomad_simulations/schema_packages/particles_state.py +++ b/src/nomad_simulations/schema_packages/particles_state.py @@ -144,79 +144,79 @@ def __init__( particles = None - # if hasattr(symbols, 'get_positions'): - # atoms = symbols - # symbols = None - # elif ( - # isinstance(symbols, (list, tuple)) - # and len(symbols) > 0 - # and isinstance(symbols[0], Atom) - # ): - # # Get data from a list or tuple of Atom objects: - # data = [ - # [atom.get_raw(name) for atom in symbols] - # for name in [ - # 'position', - # 'number', - # 'tag', - # 'momentum', - # 'mass', - # 'magmom', - # 'charge', - # ] - # ] - # atoms = self.__class__(None, *data) - # symbols = None - - # if atoms is not None: - # # Get data from another Atoms object: - # if scaled_positions is not None: - # raise NotImplementedError - # if symbols is None and numbers is None: - # numbers = atoms.get_atomic_numbers() - # if positions is None: - # positions = atoms.get_positions() - # if tags is None and atoms.has('tags'): - # tags = atoms.get_tags() - # if momenta is None and atoms.has('momenta'): - # momenta = atoms.get_momenta() - # if magmoms is None and atoms.has('initial_magmoms'): - # magmoms = atoms.get_initial_magnetic_moments() - # if masses is None and atoms.has('masses'): - # masses = atoms.get_masses() - # if charges is None and atoms.has('initial_charges'): - # charges = atoms.get_initial_charges() - # if cell is None: - # cell = atoms.get_cell() - # if celldisp is None: - # celldisp = atoms.get_celldisp() - # if pbc is None: - # pbc = atoms.get_pbc() - - # self.arrays = {} - - # if symbols is None: - # if numbers is None: - # if positions is not None: - # natoms = len(positions) - # elif scaled_positions is not None: - # natoms = len(scaled_positions) - # else: - # natoms = 0 - # numbers = np.zeros(natoms, int) - # self.new_array('numbers', numbers, int) - # else: - # if numbers is not None: - # raise TypeError('Use only one of "symbols" and "numbers".') - # else: - # self.new_array('numbers', symbols2numbers(symbols), int) - - # if self.numbers.ndim != 1: - # raise ValueError('"numbers" must be 1-dimensional.') - - # if cell is None: - # cell = np.zeros((3, 3)) - # self.set_cell(cell) + # if hasattr(types, 'get_positions'): + # atoms = types + # types = None + # elif ( + # isinstance(types, (list, tuple)) + # and len(types) > 0 + # and isinstance(types[0], Atom) + # ): + # # Get data from a list or tuple of Atom objects: + # data = [ + # [atom.get_raw(name) for atom in types] + # for name in [ + # 'position', + # 'number', + # 'tag', + # 'momentum', + # 'mass', + # 'magmom', + # 'charge', + # ] + # ] + # atoms = self.__class__(None, *data) + # types = None + + # if atoms is not None: + # # Get data from another Atoms object: + # if scaled_positions is not None: + # raise NotImplementedError + # if types is None and numbers is None: + # numbers = atoms.get_atomic_numbers() + # if positions is None: + # positions = atoms.get_positions() + # if tags is None and atoms.has('tags'): + # tags = atoms.get_tags() + # if momenta is None and atoms.has('momenta'): + # momenta = atoms.get_momenta() + # if magmoms is None and atoms.has('initial_magmoms'): + # magmoms = atoms.get_initial_magnetic_moments() + # if masses is None and atoms.has('masses'): + # masses = atoms.get_masses() + # if charges is None and atoms.has('initial_charges'): + # charges = atoms.get_initial_charges() + # if cell is None: + # cell = atoms.get_cell() + # if celldisp is None: + # celldisp = atoms.get_celldisp() + # if pbc is None: + # pbc = atoms.get_pbc() + + # self.arrays = {} + + # if types is None: + # if numbers is None: + # if positions is not None: + # natoms = len(positions) + # elif scaled_positions is not None: + # natoms = len(scaled_positions) + # else: + # natoms = 0 + # numbers = np.zeros(natoms, int) + # self.new_array('numbers', numbers, int) + # else: + # if numbers is not None: + # raise TypeError('Use only one of "types" and "numbers".') + # else: + # self.new_array('numbers', types2numbers(types), int) + + # if self.numbers.ndim != 1: + # raise ValueError('"numbers" must be 1-dimensional.') + + if cell is None: + cell = np.zeros((3, 3)) + self.set_cell(cell) # if celldisp is None: # celldisp = np.zeros(shape=(3, 1)) @@ -230,7 +230,7 @@ def __init__( # positions = np.dot(scaled_positions, self.cell) # else: # if scaled_positions is not None: - # raise TypeError('Use only one of "symbols" and "numbers".') + # raise TypeError('Use only one of "types" and "numbers".') # self.new_array('positions', positions, float, (3,)) # self.set_tags(default(tags, 0)) # self.set_masses(default(masses, None)) @@ -336,152 +336,6 @@ def get_pbc(self): """Get periodic boundary condition flags.""" return self.pbc.copy() - def new_array(self, name, a, dtype=None, shape=None): - """Add new array. - - If *shape* is not *None*, the shape of *a* will be checked.""" - - if dtype is not None: - a = np.array(a, dtype, order='C') - if len(a) == 0 and shape is not None: - a.shape = (-1,) + shape - else: - if not a.flags['C_CONTIGUOUS']: - a = np.ascontiguousarray(a) - else: - a = a.copy() - - if name in self.arrays: - raise RuntimeError(f'Array {name} already present') - - for b in self.arrays.values(): - if len(a) != len(b): - raise ValueError( - 'Array "%s" has wrong length: %d != %d.' % (name, len(a), len(b)) - ) - break - - if shape is not None and a.shape[1:] != shape: - raise ValueError( - 'Array "%s" has wrong shape %s != %s.' - % (name, a.shape, (a.shape[0:1] + shape)) - ) - - self.arrays[name] = a - - def get_array(self, name, copy=True): - """Get an array. - - Returns a copy unless the optional argument copy is false. - """ - if copy: - return self.arrays[name].copy() - else: - return self.arrays[name] - - def set_array(self, name, a, dtype=None, shape=None): - """Update array. - - If *shape* is not *None*, the shape of *a* will be checked. - If *a* is *None*, then the array is deleted.""" - - b = self.arrays.get(name) - if b is None: - if a is not None: - self.new_array(name, a, dtype, shape) - else: - if a is None: - del self.arrays[name] - else: - a = np.asarray(a) - if a.shape != b.shape: - raise ValueError( - 'Array "%s" has wrong shape %s != %s.' - % (name, a.shape, b.shape) - ) - b[:] = a - - def has(self, name): - """Check for existence of array. - - name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms', - 'initial_charges'.""" - # XXX extend has to calculator properties - return name in self.arrays - - def set_tags(self, tags): - """Set tags for all atoms. If only one tag is supplied, it is - applied to all atoms.""" - if isinstance(tags, int): - tags = [tags] * len(self) - self.set_array('tags', tags, int, ()) - - def get_tags(self): - """Get integer array of tags.""" - if 'tags' in self.arrays: - return self.arrays['tags'].copy() - else: - return np.zeros(len(self), int) - - def set_momenta(self, momenta, apply_constraint=True): - """Set momenta.""" - if apply_constraint and len(self.constraints) > 0 and momenta is not None: - momenta = np.array(momenta) # modify a copy - for constraint in self.constraints: - if hasattr(constraint, 'adjust_momenta'): - constraint.adjust_momenta(self, momenta) - self.set_array('momenta', momenta, float, (3,)) - - def set_velocities(self, velocities): - """Set the momenta by specifying the velocities.""" - self.set_momenta(self.get_masses()[:, np.newaxis] * velocities) - - def get_momenta(self): - """Get array of momenta.""" - if 'momenta' in self.arrays: - return self.arrays['momenta'].copy() - else: - return np.zeros((len(self), 3)) - - # def set_masses(self, masses='defaults'): - # """Set atomic masses in atomic mass units. - - # The array masses should contain a list of masses. In case - # the masses argument is not given or for those elements of the - # masses list that are None, standard values are set.""" - - # if isinstance(masses, str): - # if masses == 'defaults': - # masses = atomic_masses[self.arrays['numbers']] - # elif masses == 'most_common': - # masses = atomic_masses_common[self.arrays['numbers']] - # elif masses is None: - # pass - # elif not isinstance(masses, np.ndarray): - # masses = list(masses) - # for i, mass in enumerate(masses): - # if mass is None: - # masses[i] = atomic_masses[self.numbers[i]] - # self.set_array('masses', masses, float, ()) - - # def get_masses(self): - # """Get array of masses in atomic mass units.""" - # if 'masses' in self.arrays: - # return self.arrays['masses'].copy() - # else: - # return atomic_masses[self.arrays['numbers']] - - def get_charges(self): - """Get calculated charges.""" - if self._calc is None: - raise RuntimeError('Atoms object has no calculator.') - try: - return self._calc.get_charges(self) - except AttributeError: - from ase.calculators.calculator import PropertyNotImplementedError - - raise PropertyNotImplementedError - def set_positions(self, newpositions, apply_constraint=True): """Set positions, honoring any constraints. To ignore constraints, use *apply_constraint=False*.""" @@ -554,17 +408,6 @@ def wrap(self, **wrap_kw): self.positions[:] = self.get_positions(wrap=True, **wrap_kw) - # @deprecated('Please use atoms.cell.volume') - # We kind of want to deprecate this, but the ValueError behaviour - # might be desirable. Should we do this? - def get_volume(self): - """Get volume of unit cell.""" - if self.cell.rank != 3: - raise ValueError( - f'You have {self.cell.rank} lattice vectors: volume not defined' - ) - return self.cell.volume - def _get_positions(self): """Return reference to positions-array for in-place manipulations.""" return self.arrays['positions'] From c84d77e89f75110dcc4d1b6a7b2976f7a36a0344 Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Wed, 27 Nov 2024 16:08:21 +0100 Subject: [PATCH 15/16] fix rebase --- .../schema_packages/model_system.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index 7f47f0bf..f23c3e19 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -15,8 +15,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # - import re +import sys from functools import lru_cache from hashlib import sha1 from typing import TYPE_CHECKING, Optional @@ -51,7 +51,7 @@ from structlog.stdlib import BoundLogger from nomad_simulations.schema_packages.atoms_state import AtomsState -from nomad_simulations.schema_packages.particles_state import Particles, Particles, ParticlesState +from nomad_simulations.schema_packages.particles_state import Particles, ParticlesState from nomad_simulations.schema_packages.utils import ( catch_not_implemented, get_sibling_section, @@ -660,11 +660,11 @@ def to_particles(self, logger: 'BoundLogger') -> Optional[Particles]: self.periodic_boundary_conditions = [False, False, False] particles.set_pbc(pbc=self.periodic_boundary_conditions) - # # Lattice vectors - # if self.lattice_vectors is not None: - # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) - # else: - # logger.info('Could not find `AtomicCell.lattice_vectors`.') + # # Lattice vectors + # if self.lattice_vectors is not None: + # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + # else: + # logger.info('Could not find `AtomicCell.lattice_vectors`.') # Positions if self.positions is not None: @@ -1211,7 +1211,6 @@ class ModelSystem(System): ) cell = SubSection(sub_section=Cell.m_def, repeats=True) - print(f'cell: {Cell.m_def}') symmetry = SubSection(sub_section=Symmetry.m_def, repeats=True) From 243cbb6d4454305be00e08b86b09df61a28cc7d6 Mon Sep 17 00:00:00 2001 From: Bernadette Mohr Date: Tue, 3 Dec 2024 17:29:49 +0100 Subject: [PATCH 16/16] marked issues from calls to ase functions --- .../schema_packages/general.py | 2 ++ .../schema_packages/model_system.py | 30 +++++++++++-------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/src/nomad_simulations/schema_packages/general.py b/src/nomad_simulations/schema_packages/general.py index 9a2d48f0..59777858 100644 --- a/src/nomad_simulations/schema_packages/general.py +++ b/src/nomad_simulations/schema_packages/general.py @@ -217,6 +217,7 @@ def _set_system_branch_depth( system_parent=system_child, branch_depth=branch_depth + 1 ) + #! Generalize from checks for atomic systems, error with CG input def resolve_composition_formula(self, system_parent: ModelSystem) -> None: """Determine and set the composition formula for `system_parent` and all of its descendants. @@ -275,6 +276,7 @@ def get_composition_recurs(system: ModelSystem, atom_labels: list[str]) -> None: for subsystem in subsystems: get_composition_recurs(system=subsystem, atom_labels=atom_labels) + # ! CG: system_parent.cell[0].particles_state instead of atoms_state! atoms_state = ( system_parent.cell[0].atoms_state if system_parent.cell is not None else [] ) diff --git a/src/nomad_simulations/schema_packages/model_system.py b/src/nomad_simulations/schema_packages/model_system.py index f23c3e19..21cd6ad4 100644 --- a/src/nomad_simulations/schema_packages/model_system.py +++ b/src/nomad_simulations/schema_packages/model_system.py @@ -305,6 +305,7 @@ class Cell(GeometricSpace): """, ) + # TODO: default "unavailable"? type = Quantity( type=MEnum('original', 'primitive', 'conventional'), description=""" @@ -660,11 +661,11 @@ def to_particles(self, logger: 'BoundLogger') -> Optional[Particles]: self.periodic_boundary_conditions = [False, False, False] particles.set_pbc(pbc=self.periodic_boundary_conditions) - # # Lattice vectors - # if self.lattice_vectors is not None: - # ase_atoms.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) - # else: - # logger.info('Could not find `AtomicCell.lattice_vectors`.') + # Lattice vectors + if self.lattice_vectors is not None: + particles.set_cell(cell=self.lattice_vectors.to('angstrom').magnitude) + else: + logger.info('Could not find `ParticleCell.lattice_vectors`.') # Positions if self.positions is not None: @@ -1088,6 +1089,7 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: self.m_cache['elemental_composition'] = formula.elemental_composition() +# TODO: generalize! indices instead of atom_indices, state instead of atoms_state... class ModelSystem(System): """ Model system used as an input for simulating the material. @@ -1252,6 +1254,7 @@ class ModelSystem(System): """, ) + # TODO: make this work with non_atomic systems: global_composition_formula of entire system with respect to lower layers composition_formula = Quantity( type=str, description=""" @@ -1367,11 +1370,12 @@ def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None: sec_symmetry = self.m_create(Symmetry) sec_symmetry.normalize(archive, logger) - # Creating and normalizing ChemicalFormula section - # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) - sec_chemical_formula = self.m_create(ChemicalFormula) - sec_chemical_formula.normalize(archive, logger) - if sec_chemical_formula.m_cache: - self.elemental_composition = sec_chemical_formula.m_cache.get( - 'elemental_composition', [] - ) + #! ChemicalFormula calls `ase_atoms = atomic_cell.to_ase_atoms(logger=logger)` and `ase_atoms.get_chemical_formula()` + # Creating and normalizing ChemicalFormula section + # TODO add support for fractional formulas (possibly add `AtomicCell.concentrations` for each species) + sec_chemical_formula = self.m_create(ChemicalFormula) + sec_chemical_formula.normalize(archive, logger) + if sec_chemical_formula.m_cache: + self.elemental_composition = sec_chemical_formula.m_cache.get( + 'elemental_composition', [] + )