Coverage for /builds/ase/ase/ase/atoms.py: 93.71%
985 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3# Copyright 2008, 2009 CAMd
4# (see accompanying license files for details).
6"""Definition of the Atoms class.
8This module defines the central object in the ASE package: the Atoms
9object.
10"""
11from __future__ import annotations
13import copy
14import numbers
15from math import cos, pi, sin
16from typing import Sequence, Union, overload
18import numpy as np
20from ase import units
21from ase.atom import Atom
22from ase.cell import Cell
23from ase.data import atomic_masses, atomic_masses_common
24from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
25from ase.symbols import Symbols, symbols2numbers
26from ase.utils import deprecated, string2index
29class Atoms:
30 """Atoms object.
32 The Atoms object can represent an isolated molecule, or a
33 periodically repeated structure. It has a unit cell and
34 there may be periodic boundary conditions along any of the three
35 unit cell axes.
36 Information about the atoms (atomic numbers and position) is
37 stored in ndarrays. Optionally, there can be information about
38 tags, momenta, masses, magnetic moments and charges.
40 In order to calculate energies, forces and stresses, a calculator
41 object has to attached to the atoms object.
43 Parameters
44 ----------
45 symbols : str | list[str] | list[Atom]
46 Chemical formula, a list of chemical symbols, or list of
47 :class:`~ase.Atom` objects (mutually exclusive with ``numbers``).
49 - ``'H2O'``
50 - ``'COPt12'``
51 - ``['H', 'H', 'O']``
52 - ``[Atom('Ne', (x, y, z)), ...]``
54 positions : list[tuple[float, float, float]]
55 Atomic positions in Cartesian coordinates
56 (mutually exclusive with ``scaled_positions``).
57 Anything that can be converted to an ndarray of shape (n, 3) works:
58 [(x0, y0, z0), (x1, y1, z1), ...].
59 scaled_positions : list[tuple[float, float, float]]
60 Atomic positions in units of the unit cell
61 (mutually exclusive with ``positions``).
62 numbers : list[int]
63 Atomic numbers (mutually exclusive with ``symbols``).
64 tags : list[int]
65 Special purpose tags.
66 momenta : list[tuple[float, float, float]]
67 Momenta for all atoms in Cartesian coordinates
68 (mutually exclusive with ``velocities``).
69 velocities : list[tuple[float, float, float]]
70 Velocities for all atoms in Cartesian coordinates
71 (mutually exclusive with ``momenta``).
72 masses : list[float]
73 Atomic masses in atomic units.
74 magmoms : list[float] | list[tuple[float, float, float]]
75 Magnetic moments. Can be either a single value for each atom
76 for collinear calculations or three numbers for each atom for
77 non-collinear calculations.
78 charges : list[float]
79 Initial atomic charges.
80 cell : 3x3 matrix or length 3 or 6 vector, default: (0, 0, 0)
81 Unit cell vectors. Can also be given as just three
82 numbers for orthorhombic cells, or 6 numbers, where
83 first three are lengths of unit cell vectors, and the
84 other three are angles between them (in degrees), in following order:
85 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)].
86 First vector will lie in x-direction, second in xy-plane,
87 and the third one in z-positive subspace.
88 celldisp : tuple[float, float, float], default: (0, 0, 0)
89 Unit cell displacement vector. To visualize a displaced cell
90 around the center of mass of a Systems of atoms.
91 pbc : bool | tuple[bool, bool, bool], default: False
92 Periodic boundary conditions flags.
94 - ``True``
95 - ``False``
96 - ``0``
97 - ``1``
98 - ``(1, 1, 0)``
99 - ``(True, False, False)``
101 constraint : constraint object(s)
102 One or more ASE constraints applied during structure optimization.
103 calculator : calculator object
104 ASE calculator to obtain energies and atomic forces.
105 info : dict | None, default: None
106 Dictionary with additional information about the system.
107 The following keys may be used by ASE:
109 - spacegroup: :class:`~ase.spacegroup.Spacegroup` instance
110 - unit_cell: 'conventional' | 'primitive' | int | 3 ints
111 - adsorbate_info: Information about special adsorption sites
113 Items in the info attribute survives copy and slicing and can
114 be stored in and retrieved from trajectory files given that the
115 key is a string, the value is JSON-compatible and, if the value is a
116 user-defined object, its base class is importable. One should
117 not make any assumptions about the existence of keys.
119 Examples
120 --------
121 >>> from ase import Atom
123 N2 molecule (These three are equivalent):
125 >>> d = 1.104 # N2 bondlength
126 >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)])
127 >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)])
128 >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))])
130 FCC gold:
132 >>> a = 4.05 # Gold lattice constant
133 >>> b = a / 2
134 >>> fcc = Atoms('Au',
135 ... cell=[(0, b, b), (b, 0, b), (b, b, 0)],
136 ... pbc=True)
138 Hydrogen wire:
140 >>> d = 0.9 # H-H distance
141 >>> h = Atoms('H', positions=[(0, 0, 0)],
142 ... cell=(d, 0, 0),
143 ... pbc=(1, 0, 0))
144 """
146 ase_objtype = 'atoms' # For JSONability
148 def __init__(self, symbols=None,
149 positions=None, numbers=None,
150 tags=None, momenta=None, masses=None,
151 magmoms=None, charges=None,
152 scaled_positions=None,
153 cell=None, pbc=None, celldisp=None,
154 constraint=None,
155 calculator=None,
156 info=None,
157 velocities=None):
159 self._cellobj = Cell.new()
160 self._pbc = np.zeros(3, bool)
162 atoms = None
164 if hasattr(symbols, 'get_positions'):
165 atoms = symbols
166 symbols = None
167 elif (isinstance(symbols, (list, tuple)) and
168 len(symbols) > 0 and isinstance(symbols[0], Atom)):
169 # Get data from a list or tuple of Atom objects:
170 data = [[atom.get_raw(name) for atom in symbols]
171 for name in
172 ['position', 'number', 'tag', 'momentum',
173 'mass', 'magmom', 'charge']]
174 atoms = self.__class__(None, *data)
175 symbols = None
177 if atoms is not None:
178 # Get data from another Atoms object:
179 if scaled_positions is not None:
180 raise NotImplementedError
181 if symbols is None and numbers is None:
182 numbers = atoms.get_atomic_numbers()
183 if positions is None:
184 positions = atoms.get_positions()
185 if tags is None and atoms.has('tags'):
186 tags = atoms.get_tags()
187 if momenta is None and atoms.has('momenta'):
188 momenta = atoms.get_momenta()
189 if magmoms is None and atoms.has('initial_magmoms'):
190 magmoms = atoms.get_initial_magnetic_moments()
191 if masses is None and atoms.has('masses'):
192 masses = atoms.get_masses()
193 if charges is None and atoms.has('initial_charges'):
194 charges = atoms.get_initial_charges()
195 if cell is None:
196 cell = atoms.get_cell()
197 if celldisp is None:
198 celldisp = atoms.get_celldisp()
199 if pbc is None:
200 pbc = atoms.get_pbc()
201 if constraint is None:
202 constraint = [c.copy() for c in atoms.constraints]
203 if calculator is None:
204 calculator = atoms.calc
205 if info is None:
206 info = copy.deepcopy(atoms.info)
208 self.arrays: dict[str, np.ndarray] = {}
210 if symbols is not None and numbers is not None:
211 raise TypeError('Use only one of "symbols" and "numbers".')
212 if symbols is not None:
213 numbers = symbols2numbers(symbols)
214 elif numbers is None:
215 if positions is not None:
216 natoms = len(positions)
217 elif scaled_positions is not None:
218 natoms = len(scaled_positions)
219 else:
220 natoms = 0
221 numbers = np.zeros(natoms, int)
222 self.new_array('numbers', numbers, int)
224 if self.numbers.ndim != 1:
225 raise ValueError('"numbers" must be 1-dimensional.')
227 if cell is None:
228 cell = np.zeros((3, 3))
229 self.set_cell(cell)
231 if celldisp is None:
232 celldisp = np.zeros(shape=(3, 1))
233 self.set_celldisp(celldisp)
235 if positions is None:
236 if scaled_positions is None:
237 positions = np.zeros((len(self.arrays['numbers']), 3))
238 else:
239 assert self.cell.rank == 3
240 positions = np.dot(scaled_positions, self.cell)
241 else:
242 if scaled_positions is not None:
243 raise TypeError(
244 'Use only one of "positions" and "scaled_positions".')
245 self.new_array('positions', positions, float, (3,))
247 self.set_constraint(constraint)
248 self.set_tags(default(tags, 0))
249 self.set_masses(default(masses, None))
250 self.set_initial_magnetic_moments(default(magmoms, 0.0))
251 self.set_initial_charges(default(charges, 0.0))
252 if pbc is None:
253 pbc = False
254 self.set_pbc(pbc)
255 self.set_momenta(default(momenta, (0.0, 0.0, 0.0)),
256 apply_constraint=False)
258 if velocities is not None:
259 if momenta is None:
260 self.set_velocities(velocities)
261 else:
262 raise TypeError(
263 'Use only one of "momenta" and "velocities".')
265 if info is None:
266 self.info = {}
267 else:
268 self.info = dict(info)
270 self.calc = calculator
272 @property
273 def symbols(self):
274 """Get chemical symbols as a :class:`ase.symbols.Symbols` object.
276 The object works like ``atoms.numbers`` except its values
277 are strings. It supports in-place editing."""
278 return Symbols(self.numbers)
280 @symbols.setter
281 def symbols(self, obj):
282 new_symbols = Symbols.fromsymbols(obj)
283 self.numbers[:] = new_symbols.numbers
285 def get_chemical_symbols(self):
286 """Get list of chemical symbol strings.
288 Equivalent to ``list(atoms.symbols)``."""
289 return list(self.symbols)
291 def set_chemical_symbols(self, symbols):
292 """Set chemical symbols."""
293 self.set_array('numbers', symbols2numbers(symbols), int, ())
295 @property
296 def numbers(self):
297 """Attribute for direct manipulation of the atomic numbers."""
298 return self.arrays['numbers']
300 @numbers.setter
301 def numbers(self, numbers):
302 self.set_atomic_numbers(numbers)
304 def set_atomic_numbers(self, numbers):
305 """Set atomic numbers."""
306 self.set_array('numbers', numbers, int, ())
308 def get_atomic_numbers(self):
309 """Get integer array of atomic numbers."""
310 return self.arrays['numbers'].copy()
312 @property
313 def positions(self):
314 """Attribute for direct manipulation of the positions."""
315 return self.arrays['positions']
317 @positions.setter
318 def positions(self, pos):
319 self.arrays['positions'][:] = pos
321 def set_positions(self, newpositions, apply_constraint=True):
322 """Set positions, honoring any constraints. To ignore constraints,
323 use *apply_constraint=False*."""
324 if self.constraints and apply_constraint:
325 newpositions = np.array(newpositions, float)
326 for constraint in self.constraints:
327 constraint.adjust_positions(self, newpositions)
329 self.set_array('positions', newpositions, shape=(3,))
331 def get_positions(self, wrap=False, **wrap_kw):
332 """Get array of positions.
334 Parameters:
336 wrap: bool
337 wrap atoms back to the cell before returning positions
338 wrap_kw: (keyword=value) pairs
339 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
340 see :func:`ase.geometry.wrap_positions`
341 """
342 from ase.geometry import wrap_positions
343 if wrap:
344 if 'pbc' not in wrap_kw:
345 wrap_kw['pbc'] = self.pbc
346 return wrap_positions(self.positions, self.cell, **wrap_kw)
347 else:
348 return self.arrays['positions'].copy()
350 @deprecated("Please use atoms.calc = calc", FutureWarning)
351 def set_calculator(self, calc=None):
352 """Attach calculator object.
354 .. deprecated:: 3.20.0
355 Please use the equivalent ``atoms.calc = calc`` instead of this
356 method.
357 """
359 self.calc = calc
361 @deprecated("Please use atoms.calc", FutureWarning)
362 def get_calculator(self):
363 """Get currently attached calculator object.
365 .. deprecated:: 3.20.0
366 Please use the equivalent ``atoms.calc`` instead of
367 ``atoms.get_calculator()``.
368 """
370 return self.calc
372 @property
373 def calc(self):
374 """Calculator object."""
375 return self._calc
377 @calc.setter
378 def calc(self, calc):
379 self._calc = calc
380 if hasattr(calc, 'set_atoms'):
381 calc.set_atoms(self)
383 @calc.deleter
384 @deprecated('Please use atoms.calc = None', FutureWarning)
385 def calc(self):
386 """Delete calculator
388 .. deprecated:: 3.20.0
389 Please use ``atoms.calc = None``
390 """
391 self._calc = None
393 @property
394 @deprecated('Please use atoms.cell.rank instead', DeprecationWarning)
395 def number_of_lattice_vectors(self):
396 """Number of (non-zero) lattice vectors.
398 .. deprecated:: 3.21.0
399 Please use ``atoms.cell.rank`` instead
400 """
401 return self.cell.rank
403 @property
404 def constraints(self):
405 return self._constraints
407 @constraints.setter
408 def constraints(self, constraint):
409 self.set_constraint(constraint)
411 @constraints.deleter
412 def constraints(self):
413 self._constraints = []
415 def set_constraint(self, constraint=None):
416 """Apply one or more constrains.
418 The *constraint* argument must be one constraint object or a
419 list of constraint objects."""
420 if constraint is None:
421 self._constraints = []
422 else:
423 if isinstance(constraint, list):
424 self._constraints = constraint
425 elif isinstance(constraint, tuple):
426 self._constraints = list(constraint)
427 else:
428 self._constraints = [constraint]
430 def get_number_of_degrees_of_freedom(self):
431 """Calculate the number of degrees of freedom in the system."""
432 return len(self) * 3 - sum(
433 c.get_removed_dof(self) for c in self._constraints
434 )
436 def set_cell(self, cell, scale_atoms=False, apply_constraint=True):
437 """Set unit cell vectors.
439 Parameters:
441 cell: 3x3 matrix or length 3 or 6 vector
442 Unit cell. A 3x3 matrix (the three unit cell vectors) or
443 just three numbers for an orthorhombic cell. Another option is
444 6 numbers, which describes unit cell with lengths of unit cell
445 vectors and with angles between them (in degrees), in following
446 order: [len(a), len(b), len(c), angle(b,c), angle(a,c),
447 angle(a,b)]. First vector will lie in x-direction, second in
448 xy-plane, and the third one in z-positive subspace.
449 scale_atoms: bool
450 Fix atomic positions or move atoms with the unit cell?
451 Default behavior is to *not* move the atoms (scale_atoms=False).
452 apply_constraint: bool
453 Whether to apply constraints to the given cell.
455 Examples:
457 Two equivalent ways to define an orthorhombic cell:
459 >>> atoms = Atoms('He')
460 >>> a, b, c = 7, 7.5, 8
461 >>> atoms.set_cell([a, b, c])
462 >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)])
464 FCC unit cell:
466 >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)])
468 Hexagonal unit cell:
470 >>> atoms.set_cell([a, a, c, 90, 90, 120])
472 Rhombohedral unit cell:
474 >>> alpha = 77
475 >>> atoms.set_cell([a, a, a, alpha, alpha, alpha])
476 """
478 # Override pbcs if and only if given a Cell object:
479 cell = Cell.new(cell)
481 # XXX not working well during initialize due to missing _constraints
482 if apply_constraint and hasattr(self, '_constraints'):
483 for constraint in self.constraints:
484 if hasattr(constraint, 'adjust_cell'):
485 constraint.adjust_cell(self, cell)
487 if scale_atoms:
488 M = np.linalg.solve(self.cell.complete(), cell.complete())
489 self.positions[:] = np.dot(self.positions, M)
491 self.cell[:] = cell
493 def set_celldisp(self, celldisp):
494 """Set the unit cell displacement vectors."""
495 celldisp = np.array(celldisp, float)
496 self._celldisp = celldisp
498 def get_celldisp(self):
499 """Get the unit cell displacement vectors."""
500 return self._celldisp.copy()
502 def get_cell(self, complete=False):
503 """Get the three unit cell vectors as a `class`:ase.cell.Cell` object.
505 The Cell object resembles a 3x3 ndarray, and cell[i, j]
506 is the jth Cartesian coordinate of the ith cell vector."""
507 if complete:
508 cell = self.cell.complete()
509 else:
510 cell = self.cell.copy()
512 return cell
514 @deprecated('Please use atoms.cell.cellpar() instead', DeprecationWarning)
515 def get_cell_lengths_and_angles(self):
516 """Get unit cell parameters. Sequence of 6 numbers.
518 First three are unit cell vector lengths and second three
519 are angles between them::
521 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
523 in degrees.
525 .. deprecated:: 3.21.0
526 Please use ``atoms.cell.cellpar()`` instead
527 """
528 return self.cell.cellpar()
530 @deprecated('Please use atoms.cell.reciprocal()', DeprecationWarning)
531 def get_reciprocal_cell(self):
532 """Get the three reciprocal lattice vectors as a 3x3 ndarray.
534 Note that the commonly used factor of 2 pi for Fourier
535 transforms is not included here.
537 .. deprecated:: 3.21.0
538 Please use ``atoms.cell.reciprocal()``
539 """
540 return self.cell.reciprocal()
542 @property
543 def pbc(self):
544 """Reference to pbc-flags for in-place manipulations."""
545 return self._pbc
547 @pbc.setter
548 def pbc(self, pbc):
549 self._pbc[:] = pbc
551 def set_pbc(self, pbc):
552 """Set periodic boundary condition flags."""
553 self.pbc = pbc
555 def get_pbc(self):
556 """Get periodic boundary condition flags."""
557 return self.pbc.copy()
559 def new_array(self, name, a, dtype=None, shape=None):
560 """Add new array.
562 If *shape* is not *None*, the shape of *a* will be checked."""
564 if dtype is not None:
565 a = np.array(a, dtype, order='C')
566 if len(a) == 0 and shape is not None:
567 a.shape = (-1,) + shape
568 else:
569 if not a.flags['C_CONTIGUOUS']:
570 a = np.ascontiguousarray(a)
571 else:
572 a = a.copy()
574 if name in self.arrays:
575 raise RuntimeError(f'Array {name} already present')
577 for b in self.arrays.values():
578 if len(a) != len(b):
579 raise ValueError('Array "%s" has wrong length: %d != %d.' %
580 (name, len(a), len(b)))
581 break
583 if shape is not None and a.shape[1:] != shape:
584 raise ValueError(
585 f'Array "{name}" has wrong shape {a.shape} != '
586 f'{(a.shape[0:1] + shape)}.')
588 self.arrays[name] = a
590 def get_array(self, name, copy=True):
591 """Get an array.
593 Returns a copy unless the optional argument copy is false.
594 """
595 if copy:
596 return self.arrays[name].copy()
597 else:
598 return self.arrays[name]
600 def set_array(self, name, a, dtype=None, shape=None):
601 """Update array.
603 If *shape* is not *None*, the shape of *a* will be checked.
604 If *a* is *None*, then the array is deleted."""
606 b = self.arrays.get(name)
607 if b is None:
608 if a is not None:
609 self.new_array(name, a, dtype, shape)
610 else:
611 if a is None:
612 del self.arrays[name]
613 else:
614 a = np.asarray(a)
615 if a.shape != b.shape:
616 raise ValueError(
617 f'Array "{name}" has wrong shape '
618 f'{a.shape} != {b.shape}.')
619 b[:] = a
621 def has(self, name):
622 """Check for existence of array.
624 name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms',
625 'initial_charges'."""
626 # XXX extend has to calculator properties
627 return name in self.arrays
629 def get_chemical_formula(self, mode='hill', empirical=False):
630 """Get the chemical formula as a string based on the chemical symbols.
632 Parameters:
634 mode: str
635 There are four different modes available:
637 'all': The list of chemical symbols are contracted to a string,
638 e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'.
640 'reduce': The same as 'all' where repeated elements are contracted
641 to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to
642 'CH3OCH3'.
644 'hill': The list of chemical symbols are contracted to a string
645 following the Hill notation (alphabetical order with C and H
646 first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to
647 'H2O4S'. This is default.
649 'metal': The list of chemical symbols (alphabetical metals,
650 and alphabetical non-metals)
652 empirical, bool (optional, default=False)
653 Divide the symbol counts by their greatest common divisor to yield
654 an empirical formula. Only for mode `metal` and `hill`.
655 """
656 return self.symbols.get_chemical_formula(mode, empirical)
658 def set_tags(self, tags):
659 """Set tags for all atoms. If only one tag is supplied, it is
660 applied to all atoms."""
661 if isinstance(tags, int):
662 tags = [tags] * len(self)
663 self.set_array('tags', tags, int, ())
665 def get_tags(self):
666 """Get integer array of tags."""
667 if 'tags' in self.arrays:
668 return self.arrays['tags'].copy()
669 else:
670 return np.zeros(len(self), int)
672 def set_momenta(self, momenta, apply_constraint=True):
673 """Set momenta."""
674 if (apply_constraint and len(self.constraints) > 0 and
675 momenta is not None):
676 momenta = np.array(momenta) # modify a copy
677 for constraint in self.constraints:
678 if hasattr(constraint, 'adjust_momenta'):
679 constraint.adjust_momenta(self, momenta)
680 self.set_array('momenta', momenta, float, (3,))
682 def get_momenta(self):
683 """Get array of momenta."""
684 if 'momenta' in self.arrays:
685 return self.arrays['momenta'].copy()
686 else:
687 return np.zeros((len(self), 3))
689 def get_velocities(self):
690 """Get array of velocities."""
691 momenta = self.get_momenta()
692 masses = self.get_masses()
693 return momenta / masses[:, np.newaxis]
695 def set_velocities(self, velocities):
696 """Set the momenta by specifying the velocities."""
697 self.set_momenta(self.get_masses()[:, np.newaxis] * velocities)
699 def set_masses(self, masses='defaults'):
700 """Set atomic masses in atomic mass units.
702 The array masses should contain a list of masses. In case
703 the masses argument is not given or for those elements of the
704 masses list that are None, standard values are set."""
706 if isinstance(masses, str):
707 if masses == 'defaults':
708 masses = atomic_masses[self.arrays['numbers']]
709 elif masses == 'most_common':
710 masses = atomic_masses_common[self.arrays['numbers']]
711 elif masses is None:
712 pass
713 elif not isinstance(masses, np.ndarray):
714 masses = list(masses)
715 for i, mass in enumerate(masses):
716 if mass is None:
717 masses[i] = atomic_masses[self.numbers[i]]
718 self.set_array('masses', masses, float, ())
720 def get_masses(self) -> np.ndarray:
721 """Get masses of atoms.
723 Returns
724 -------
725 masses : np.ndarray
726 Atomic masses in dalton (unified atomic mass units).
728 Examples
729 --------
730 >>> from ase.build import molecule
731 >>> atoms = molecule('CH4')
732 >>> atoms.get_masses()
733 array([ 12.011, 1.008, 1.008, 1.008, 1.008])
734 >>> total_mass = atoms.get_masses().sum()
735 >>> print(f'{total_mass:f}')
736 16.043000
737 """
738 if 'masses' in self.arrays:
739 return self.arrays['masses'].copy()
740 return atomic_masses[self.arrays['numbers']]
742 def set_initial_magnetic_moments(self, magmoms=None):
743 """Set the initial magnetic moments.
745 Use either one or three numbers for every atom (collinear
746 or non-collinear spins)."""
748 if magmoms is None:
749 self.set_array('initial_magmoms', None)
750 else:
751 magmoms = np.asarray(magmoms)
752 self.set_array('initial_magmoms', magmoms, float,
753 magmoms.shape[1:])
755 def get_initial_magnetic_moments(self):
756 """Get array of initial magnetic moments."""
757 if 'initial_magmoms' in self.arrays:
758 return self.arrays['initial_magmoms'].copy()
759 else:
760 return np.zeros(len(self))
762 def get_magnetic_moments(self):
763 """Get calculated local magnetic moments."""
764 if self._calc is None:
765 raise RuntimeError('Atoms object has no calculator.')
766 return self._calc.get_magnetic_moments(self)
768 def get_magnetic_moment(self):
769 """Get calculated total magnetic moment."""
770 if self._calc is None:
771 raise RuntimeError('Atoms object has no calculator.')
772 return self._calc.get_magnetic_moment(self)
774 def set_initial_charges(self, charges=None):
775 """Set the initial charges."""
777 if charges is None:
778 self.set_array('initial_charges', None)
779 else:
780 self.set_array('initial_charges', charges, float, ())
782 def get_initial_charges(self):
783 """Get array of initial charges."""
784 if 'initial_charges' in self.arrays:
785 return self.arrays['initial_charges'].copy()
786 else:
787 return np.zeros(len(self))
789 def get_charges(self):
790 """Get calculated charges."""
791 if self._calc is None:
792 raise RuntimeError('Atoms object has no calculator.')
793 try:
794 return self._calc.get_charges(self)
795 except AttributeError:
796 from ase.calculators.calculator import PropertyNotImplementedError
797 raise PropertyNotImplementedError
799 def get_potential_energy(self, force_consistent=False,
800 apply_constraint=True):
801 """Calculate potential energy.
803 Ask the attached calculator to calculate the potential energy and
804 apply constraints. Use *apply_constraint=False* to get the raw
805 forces.
807 When supported by the calculator, either the energy extrapolated
808 to zero Kelvin or the energy consistent with the forces (the free
809 energy) can be returned.
810 """
811 if self._calc is None:
812 raise RuntimeError('Atoms object has no calculator.')
813 if force_consistent:
814 energy = self._calc.get_potential_energy(
815 self, force_consistent=force_consistent)
816 else:
817 energy = self._calc.get_potential_energy(self)
818 if apply_constraint:
819 for constraint in self.constraints:
820 if hasattr(constraint, 'adjust_potential_energy'):
821 energy += constraint.adjust_potential_energy(self)
822 return energy
824 def get_properties(self, properties):
825 """This method is experimental; currently for internal use."""
826 # XXX Something about constraints.
827 if self._calc is None:
828 raise RuntimeError('Atoms object has no calculator.')
829 return self._calc.calculate_properties(self, properties)
831 def get_potential_energies(self):
832 """Calculate the potential energies of all the atoms.
834 Only available with calculators supporting per-atom energies
835 (e.g. classical potentials).
836 """
837 if self._calc is None:
838 raise RuntimeError('Atoms object has no calculator.')
839 return self._calc.get_potential_energies(self)
841 def get_kinetic_energy(self):
842 """Get the kinetic energy."""
843 momenta = self.arrays.get('momenta')
844 if momenta is None:
845 return 0.0
846 return 0.5 * np.vdot(momenta, self.get_velocities())
848 def get_total_energy(self):
849 """Get the total energy - potential plus kinetic energy."""
850 return self.get_potential_energy() + self.get_kinetic_energy()
852 def get_forces(self, apply_constraint=True, md=False):
853 """Calculate atomic forces.
855 Ask the attached calculator to calculate the forces and apply
856 constraints. Use *apply_constraint=False* to get the raw
857 forces.
859 For molecular dynamics (md=True) we don't apply the constraint
860 to the forces but to the momenta. When holonomic constraints for
861 rigid linear triatomic molecules are present, ask the constraints
862 to redistribute the forces within each triple defined in the
863 constraints (required for molecular dynamics with this type of
864 constraints)."""
866 if self._calc is None:
867 raise RuntimeError('Atoms object has no calculator.')
868 forces = self._calc.get_forces(self)
870 if apply_constraint:
871 # We need a special md flag here because for MD we want
872 # to skip real constraints but include special "constraints"
873 # Like Hookean.
874 for constraint in self.constraints:
875 if md and hasattr(constraint, 'redistribute_forces_md'):
876 constraint.redistribute_forces_md(self, forces)
877 if not md or hasattr(constraint, 'adjust_potential_energy'):
878 constraint.adjust_forces(self, forces)
879 return forces
881 # Informs calculators (e.g. Asap) that ideal gas contribution is added here.
882 _ase_handles_dynamic_stress = True
884 def get_stress(self, voigt=True, apply_constraint=True,
885 include_ideal_gas=False):
886 """Calculate stress tensor.
888 Returns an array of the six independent components of the
889 symmetric stress tensor, in the traditional Voigt order
890 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt
891 order.
893 The ideal gas contribution to the stresses is added if the
894 atoms have momenta and ``include_ideal_gas`` is set to True.
895 """
897 if self._calc is None:
898 raise RuntimeError('Atoms object has no calculator.')
900 stress = self._calc.get_stress(self)
901 shape = stress.shape
903 if shape == (3, 3):
904 # Convert to the Voigt form before possibly applying
905 # constraints and adding the dynamic part of the stress
906 # (the "ideal gas contribution").
907 stress = full_3x3_to_voigt_6_stress(stress)
908 else:
909 assert shape == (6,)
911 if apply_constraint:
912 for constraint in self.constraints:
913 if hasattr(constraint, 'adjust_stress'):
914 constraint.adjust_stress(self, stress)
916 # Add ideal gas contribution, if applicable
917 if include_ideal_gas and self.has('momenta'):
918 stress += self.get_kinetic_stress()
920 if voigt:
921 return stress
922 else:
923 return voigt_6_to_full_3x3_stress(stress)
925 def get_stresses(self, include_ideal_gas=False, voigt=True):
926 """Calculate the stress-tensor of all the atoms.
928 Only available with calculators supporting per-atom energies and
929 stresses (e.g. classical potentials). Even for such calculators
930 there is a certain arbitrariness in defining per-atom stresses.
932 The ideal gas contribution to the stresses is added if the
933 atoms have momenta and ``include_ideal_gas`` is set to True.
934 """
935 if self._calc is None:
936 raise RuntimeError('Atoms object has no calculator.')
937 stresses = self._calc.get_stresses(self)
939 # make sure `stresses` are in voigt form
940 if np.shape(stresses)[1:] == (3, 3):
941 stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses]
942 stresses = np.array(stresses_voigt)
944 # REMARK: The ideal gas contribution is intensive, i.e., the volume
945 # is divided out. We currently don't check if `stresses` are intensive
946 # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`.
947 # It might be good to check this here, but adds computational overhead.
949 if include_ideal_gas and self.has('momenta'):
950 stresses += self.get_kinetic_stresses()
952 if voigt:
953 return stresses
954 else:
955 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
956 return np.array(stresses_3x3)
958 def get_kinetic_stress(self, voigt=True):
959 """Calculate the kinetic part of the Virial stress tensor."""
960 stress = np.zeros(6) # Voigt notation
961 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
962 p = self.get_momenta()
963 masses = self.get_masses()
964 invmass = 1.0 / masses
965 invvol = 1.0 / self.get_volume()
966 for alpha in range(3):
967 for beta in range(alpha, 3):
968 stress[stresscomp[alpha, beta]] -= (
969 p[:, alpha] * p[:, beta] * invmass).sum() * invvol
971 if voigt:
972 return stress
973 else:
974 return voigt_6_to_full_3x3_stress(stress)
976 def get_kinetic_stresses(self, voigt=True):
977 """Calculate the kinetic part of the Virial stress of all the atoms."""
978 stresses = np.zeros((len(self), 6)) # Voigt notation
979 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
980 if hasattr(self._calc, 'get_atomic_volumes'):
981 invvol = 1.0 / self._calc.get_atomic_volumes()
982 else:
983 invvol = self.get_global_number_of_atoms() / self.get_volume()
984 p = self.get_momenta()
985 invmass = 1.0 / self.get_masses()
986 for alpha in range(3):
987 for beta in range(alpha, 3):
988 stresses[:, stresscomp[alpha, beta]] -= (
989 p[:, alpha] * p[:, beta] * invmass * invvol)
991 if voigt:
992 return stresses
993 else:
994 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
995 return np.array(stresses_3x3)
997 def get_dipole_moment(self):
998 """Calculate the electric dipole moment for the atoms object.
1000 Only available for calculators which has a get_dipole_moment()
1001 method."""
1003 if self._calc is None:
1004 raise RuntimeError('Atoms object has no calculator.')
1005 return self._calc.get_dipole_moment(self)
1007 def copy(self):
1008 """Return a copy."""
1009 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1010 celldisp=self._celldisp.copy())
1012 atoms.arrays = {}
1013 for name, a in self.arrays.items():
1014 atoms.arrays[name] = a.copy()
1015 atoms.constraints = copy.deepcopy(self.constraints)
1016 return atoms
1018 def todict(self):
1019 """For basic JSON (non-database) support."""
1020 d = dict(self.arrays)
1021 d['cell'] = np.asarray(self.cell)
1022 d['pbc'] = self.pbc
1023 if self._celldisp.any():
1024 d['celldisp'] = self._celldisp
1025 if self.constraints:
1026 d['constraints'] = self.constraints
1027 if self.info:
1028 d['info'] = self.info
1029 # Calculator... trouble.
1030 return d
1032 @classmethod
1033 def fromdict(cls, dct):
1034 """Rebuild atoms object from dictionary representation (todict)."""
1035 dct = dct.copy()
1036 kw = {name: dct.pop(name)
1037 for name in ['numbers', 'positions', 'cell', 'pbc']}
1038 constraints = dct.pop('constraints', None)
1039 if constraints:
1040 from ase.constraints import dict2constraint
1041 constraints = [dict2constraint(d) for d in constraints]
1043 info = dct.pop('info', None)
1045 atoms = cls(constraint=constraints,
1046 celldisp=dct.pop('celldisp', None),
1047 info=info, **kw)
1048 natoms = len(atoms)
1050 # Some arrays are named differently from the atoms __init__ keywords.
1051 # Also, there may be custom arrays. Hence we set them directly:
1052 for name, arr in dct.items():
1053 assert len(arr) == natoms, name
1054 assert isinstance(arr, np.ndarray)
1055 atoms.arrays[name] = arr
1056 return atoms
1058 def __len__(self):
1059 return len(self.arrays['positions'])
1061 @deprecated(
1062 "Please use len(self) or, if your atoms are distributed, "
1063 "self.get_global_number_of_atoms.",
1064 category=FutureWarning,
1065 )
1066 def get_number_of_atoms(self):
1067 """
1068 .. deprecated:: 3.18.1
1069 You probably want ``len(atoms)``. Or if your atoms are distributed,
1070 use (and see) :func:`get_global_number_of_atoms()`.
1071 """
1072 return len(self)
1074 def get_global_number_of_atoms(self):
1075 """Returns the global number of atoms in a distributed-atoms parallel
1076 simulation.
1078 DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING!
1080 Equivalent to len(atoms) in the standard ASE Atoms class. You should
1081 normally use len(atoms) instead. This function's only purpose is to
1082 make compatibility between ASE and Asap easier to maintain by having a
1083 few places in ASE use this function instead. It is typically only
1084 when counting the global number of degrees of freedom or in similar
1085 situations.
1086 """
1087 return len(self)
1089 def __repr__(self):
1090 tokens = []
1092 N = len(self)
1093 if N <= 60:
1094 symbols = self.get_chemical_formula('reduce')
1095 else:
1096 symbols = self.get_chemical_formula('hill')
1097 tokens.append(f"symbols='{symbols}'")
1099 if self.pbc.any() and not self.pbc.all():
1100 tokens.append(f'pbc={self.pbc.tolist()}')
1101 else:
1102 tokens.append(f'pbc={self.pbc[0]}')
1104 cell = self.cell
1105 if cell:
1106 if cell.orthorhombic:
1107 cell = cell.lengths().tolist()
1108 else:
1109 cell = cell.tolist()
1110 tokens.append(f'cell={cell}')
1112 for name in sorted(self.arrays):
1113 if name in ['numbers', 'positions']:
1114 continue
1115 tokens.append(f'{name}=...')
1117 if self.constraints:
1118 if len(self.constraints) == 1:
1119 constraint = self.constraints[0]
1120 else:
1121 constraint = self.constraints
1122 tokens.append(f'constraint={constraint!r}')
1124 if self._calc is not None:
1125 tokens.append('calculator={}(...)'
1126 .format(self._calc.__class__.__name__))
1128 return '{}({})'.format(self.__class__.__name__, ', '.join(tokens))
1130 def __add__(self, other):
1131 atoms = self.copy()
1132 atoms += other
1133 return atoms
1135 def extend(self, other):
1136 """Extend atoms object by appending atoms from *other*."""
1137 if isinstance(other, Atom):
1138 other = self.__class__([other])
1140 n1 = len(self)
1141 n2 = len(other)
1143 for name, a1 in self.arrays.items():
1144 a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype)
1145 a[:n1] = a1
1146 if name == 'masses':
1147 a2 = other.get_masses()
1148 else:
1149 a2 = other.arrays.get(name)
1150 if a2 is not None:
1151 a[n1:] = a2
1152 self.arrays[name] = a
1154 for name, a2 in other.arrays.items():
1155 if name in self.arrays:
1156 continue
1157 a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype)
1158 a[n1:] = a2
1159 if name == 'masses':
1160 a[:n1] = self.get_masses()[:n1]
1161 else:
1162 a[:n1] = 0
1164 self.set_array(name, a)
1166 def __iadd__(self, other):
1167 self.extend(other)
1168 return self
1170 def append(self, atom):
1171 """Append atom to end."""
1172 self.extend(self.__class__([atom]))
1174 def __iter__(self):
1175 for i in range(len(self)):
1176 yield self[i]
1178 @overload
1179 def __getitem__(self, i: Union[int, np.integer]) -> Atom: ...
1181 @overload
1182 def __getitem__(self, i: Union[Sequence, slice, np.ndarray]) -> Atoms: ...
1184 def __getitem__(self, i):
1185 """Return a subset of the atoms.
1187 i -- scalar integer, list of integers, or slice object
1188 describing which atoms to return.
1190 If i is a scalar, return an Atom object. If i is a list or a
1191 slice, return an Atoms object with the same cell, pbc, and
1192 other associated info as the original Atoms object. The
1193 indices of the constraints will be shuffled so that they match
1194 the indexing in the subset returned.
1196 """
1198 if isinstance(i, numbers.Integral):
1199 natoms = len(self)
1200 if i < -natoms or i >= natoms:
1201 raise IndexError('Index out of range.')
1203 return Atom(atoms=self, index=i)
1204 elif not isinstance(i, slice):
1205 i = np.array(i)
1206 if len(i) == 0:
1207 i = np.array([], dtype=int)
1208 # if i is a mask
1209 if i.dtype == bool:
1210 if len(i) != len(self):
1211 raise IndexError('Length of mask {} must equal '
1212 'number of atoms {}'
1213 .format(len(i), len(self)))
1214 i = np.arange(len(self))[i]
1216 import copy
1218 conadd = []
1219 # Constraints need to be deepcopied, but only the relevant ones.
1220 for con in copy.deepcopy(self.constraints):
1221 try:
1222 con.index_shuffle(self, i)
1223 except (IndexError, NotImplementedError):
1224 pass
1225 else:
1226 conadd.append(con)
1228 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1229 # should be communicated to the slice as well
1230 celldisp=self._celldisp)
1231 # TODO: Do we need to shuffle indices in adsorbate_info too?
1233 atoms.arrays = {}
1234 for name, a in self.arrays.items():
1235 atoms.arrays[name] = a[i].copy()
1237 atoms.constraints = conadd
1238 return atoms
1240 def __delitem__(self, i):
1241 from ase.constraints import FixAtoms
1242 for c in self._constraints:
1243 if not isinstance(c, FixAtoms):
1244 raise RuntimeError('Remove constraint using set_constraint() '
1245 'before deleting atoms.')
1247 if isinstance(i, list) and len(i) > 0:
1248 # Make sure a list of booleans will work correctly and not be
1249 # interpreted at 0 and 1 indices.
1250 i = np.array(i)
1252 if len(self._constraints) > 0:
1253 n = len(self)
1254 i = np.arange(n)[i]
1255 if isinstance(i, int):
1256 i = [i]
1257 constraints = []
1258 for c in self._constraints:
1259 c = c.delete_atoms(i, n)
1260 if c is not None:
1261 constraints.append(c)
1262 self.constraints = constraints
1264 mask = np.ones(len(self), bool)
1265 mask[i] = False
1266 for name, a in self.arrays.items():
1267 self.arrays[name] = a[mask]
1269 def pop(self, i=-1):
1270 """Remove and return atom at index *i* (default last)."""
1271 atom = self[i]
1272 atom.cut_reference_to_atoms()
1273 del self[i]
1274 return atom
1276 def __imul__(self, m):
1277 """In-place repeat of atoms."""
1278 if isinstance(m, int):
1279 m = (m, m, m)
1281 for x, vec in zip(m, self.cell):
1282 if x != 1 and not vec.any():
1283 raise ValueError('Cannot repeat along undefined lattice '
1284 'vector')
1286 M = np.prod(m)
1287 n = len(self)
1289 for name, a in self.arrays.items():
1290 self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1))
1292 positions = self.arrays['positions']
1293 i0 = 0
1294 for m0 in range(m[0]):
1295 for m1 in range(m[1]):
1296 for m2 in range(m[2]):
1297 i1 = i0 + n
1298 positions[i0:i1] += np.dot((m0, m1, m2), self.cell)
1299 i0 = i1
1301 if self.constraints is not None:
1302 self.constraints = [c.repeat(m, n) for c in self.constraints]
1304 self.cell = np.array([m[c] * self.cell[c] for c in range(3)])
1306 return self
1308 def repeat(self, rep):
1309 """Create new repeated atoms object.
1311 The *rep* argument should be a sequence of three positive
1312 integers like *(2,3,1)* or a single integer (*r*) equivalent
1313 to *(r,r,r)*."""
1315 atoms = self.copy()
1316 atoms *= rep
1317 return atoms
1319 def __mul__(self, rep):
1320 return self.repeat(rep)
1322 def translate(self, displacement):
1323 """Translate atomic positions.
1325 The displacement argument can be a float an xyz vector or an
1326 nx3 array (where n is the number of atoms)."""
1328 self.arrays['positions'] += np.array(displacement)
1330 def center(self, vacuum=None, axis=(0, 1, 2), about=None):
1331 """Center atoms in unit cell.
1333 Centers the atoms in the unit cell, so there is the same
1334 amount of vacuum on all sides.
1336 vacuum: float (default: None)
1337 If specified adjust the amount of vacuum when centering.
1338 If vacuum=10.0 there will thus be 10 Angstrom of vacuum
1339 on each side.
1340 axis: int or sequence of ints
1341 Axis or axes to act on. Default: Act on all axes.
1342 about: float or array (default: None)
1343 If specified, center the atoms about <about>.
1344 I.e., about=(0., 0., 0.) (or just "about=0.", interpreted
1345 identically), to center about the origin.
1346 """
1348 # Find the orientations of the faces of the unit cell
1349 cell = self.cell.complete()
1350 dirs = np.zeros_like(cell)
1352 lengths = cell.lengths()
1353 for i in range(3):
1354 dirs[i] = np.cross(cell[i - 1], cell[i - 2])
1355 dirs[i] /= np.linalg.norm(dirs[i])
1356 if dirs[i] @ cell[i] < 0.0:
1357 dirs[i] *= -1
1359 if isinstance(axis, int):
1360 axes = (axis,)
1361 else:
1362 axes = axis
1364 # Now, decide how much each basis vector should be made longer
1365 pos = self.positions
1366 longer = np.zeros(3)
1367 shift = np.zeros(3)
1368 for i in axes:
1369 if len(pos):
1370 scalarprod = pos @ dirs[i]
1371 p0 = scalarprod.min()
1372 p1 = scalarprod.max()
1373 else:
1374 p0 = 0
1375 p1 = 0
1376 height = cell[i] @ dirs[i]
1377 if vacuum is not None:
1378 lng = (p1 - p0 + 2 * vacuum) - height
1379 else:
1380 lng = 0.0 # Do not change unit cell size!
1381 top = lng + height - p1
1382 shf = 0.5 * (top - p0)
1383 cosphi = cell[i] @ dirs[i] / lengths[i]
1384 longer[i] = lng / cosphi
1385 shift[i] = shf / cosphi
1387 # Now, do it!
1388 translation = np.zeros(3)
1389 for i in axes:
1390 nowlen = lengths[i]
1391 if vacuum is not None:
1392 self.cell[i] = cell[i] * (1 + longer[i] / nowlen)
1393 translation += shift[i] * cell[i] / nowlen
1395 # We calculated translations using the completed cell,
1396 # so directions without cell vectors will have been centered
1397 # along a "fake" vector of length 1.
1398 # Therefore, we adjust by -0.5:
1399 if not any(self.cell[i]):
1400 translation[i] -= 0.5
1402 # Optionally, translate to center about a point in space.
1403 if about is not None:
1404 for n, vector in enumerate(self.cell):
1405 if n in axes:
1406 translation -= vector / 2.0
1407 translation[n] += about[n]
1409 self.positions += translation
1411 def get_center_of_mass(self, scaled=False, indices=None):
1412 """Get the center of mass.
1414 Parameters
1415 ----------
1416 scaled : bool
1417 If True, the center of mass in scaled coordinates is returned.
1418 indices : list | slice | str, default: None
1419 If specified, the center of mass of a subset of atoms is returned.
1420 """
1421 if indices is None:
1422 indices = slice(None)
1423 elif isinstance(indices, str):
1424 indices = string2index(indices)
1426 masses = self.get_masses()[indices]
1427 com = masses @ self.positions[indices] / masses.sum()
1428 if scaled:
1429 return self.cell.scaled_positions(com)
1430 return com # Cartesian coordinates
1432 def set_center_of_mass(self, com, scaled=False):
1433 """Set the center of mass.
1435 If scaled=True the center of mass is expected in scaled coordinates.
1436 Constraints are considered for scaled=False.
1437 """
1438 old_com = self.get_center_of_mass(scaled=scaled)
1439 difference = com - old_com
1440 if scaled:
1441 self.set_scaled_positions(self.get_scaled_positions() + difference)
1442 else:
1443 self.set_positions(self.get_positions() + difference)
1445 def get_moments_of_inertia(self, vectors=False):
1446 """Get the moments of inertia along the principal axes.
1448 The three principal moments of inertia are computed from the
1449 eigenvalues of the symmetric inertial tensor. Periodic boundary
1450 conditions are ignored. Units of the moments of inertia are
1451 amu*angstrom**2.
1452 """
1453 com = self.get_center_of_mass()
1454 positions = self.get_positions()
1455 positions -= com # translate center of mass to origin
1456 masses = self.get_masses()
1458 # Initialize elements of the inertial tensor
1459 I11 = I22 = I33 = I12 = I13 = I23 = 0.0
1460 for i in range(len(self)):
1461 x, y, z = positions[i]
1462 m = masses[i]
1464 I11 += m * (y ** 2 + z ** 2)
1465 I22 += m * (x ** 2 + z ** 2)
1466 I33 += m * (x ** 2 + y ** 2)
1467 I12 += -m * x * y
1468 I13 += -m * x * z
1469 I23 += -m * y * z
1471 Itensor = np.array([[I11, I12, I13],
1472 [I12, I22, I23],
1473 [I13, I23, I33]])
1475 evals, evecs = np.linalg.eigh(Itensor)
1476 if vectors:
1477 return evals, evecs.transpose()
1478 else:
1479 return evals
1481 def get_angular_momentum(self):
1482 """Get total angular momentum with respect to the center of mass."""
1483 com = self.get_center_of_mass()
1484 positions = self.get_positions()
1485 positions -= com # translate center of mass to origin
1486 return np.cross(positions, self.get_momenta()).sum(0)
1488 def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False):
1489 """Rotate atoms based on a vector and an angle, or two vectors.
1491 Parameters:
1493 a = None:
1494 Angle that the atoms is rotated around the vector 'v'. 'a'
1495 can also be a vector and then 'a' is rotated
1496 into 'v'.
1498 v:
1499 Vector to rotate the atoms around. Vectors can be given as
1500 strings: 'x', '-x', 'y', ... .
1502 center = (0, 0, 0):
1503 The center is kept fixed under the rotation. Use 'COM' to fix
1504 the center of mass, 'COP' to fix the center of positions or
1505 'COU' to fix the center of cell.
1507 rotate_cell = False:
1508 If true the cell is also rotated.
1510 Examples:
1512 Rotate 90 degrees around the z-axis, so that the x-axis is
1513 rotated into the y-axis:
1515 >>> atoms = Atoms()
1516 >>> atoms.rotate(90, 'z')
1517 >>> atoms.rotate(90, (0, 0, 1))
1518 >>> atoms.rotate(-90, '-z')
1519 >>> atoms.rotate('x', 'y')
1520 >>> atoms.rotate((1, 0, 0), (0, 1, 0))
1521 """
1523 if not isinstance(a, numbers.Real):
1524 a, v = v, a
1526 norm = np.linalg.norm
1527 v = string2vector(v)
1529 normv = norm(v)
1531 if normv == 0.0:
1532 raise ZeroDivisionError('Cannot rotate: norm(v) == 0')
1534 if isinstance(a, numbers.Real):
1535 a *= pi / 180
1536 v /= normv
1537 c = cos(a)
1538 s = sin(a)
1539 else:
1540 v2 = string2vector(a)
1541 v /= normv
1542 normv2 = np.linalg.norm(v2)
1543 if normv2 == 0:
1544 raise ZeroDivisionError('Cannot rotate: norm(a) == 0')
1545 v2 /= norm(v2)
1546 c = np.dot(v, v2)
1547 v = np.cross(v, v2)
1548 s = norm(v)
1549 # In case *v* and *a* are parallel, np.cross(v, v2) vanish
1550 # and can't be used as a rotation axis. However, in this
1551 # case any rotation axis perpendicular to v2 will do.
1552 eps = 1e-7
1553 if s < eps:
1554 v = np.cross((0, 0, 1), v2)
1555 if norm(v) < eps:
1556 v = np.cross((1, 0, 0), v2)
1557 assert norm(v) >= eps
1558 elif s > 0:
1559 v /= s
1561 center = self._centering_as_array(center)
1563 p = self.arrays['positions'] - center
1564 self.arrays['positions'][:] = (c * p -
1565 np.cross(p, s * v) +
1566 np.outer(np.dot(p, v), (1.0 - c) * v) +
1567 center)
1568 if rotate_cell:
1569 rotcell = self.get_cell()
1570 rotcell[:] = (c * rotcell -
1571 np.cross(rotcell, s * v) +
1572 np.outer(np.dot(rotcell, v), (1.0 - c) * v))
1573 self.set_cell(rotcell)
1575 def _centering_as_array(self, center):
1576 if isinstance(center, str):
1577 if center.lower() == 'com':
1578 center = self.get_center_of_mass()
1579 elif center.lower() == 'cop':
1580 center = self.get_positions().mean(axis=0)
1581 elif center.lower() == 'cou':
1582 center = self.get_cell().sum(axis=0) / 2
1583 else:
1584 raise ValueError('Cannot interpret center')
1585 else:
1586 center = np.array(center, float)
1587 return center
1589 def euler_rotate(
1590 self,
1591 phi: float = 0.0,
1592 theta: float = 0.0,
1593 psi: float = 0.0,
1594 center: Sequence[float] = (0.0, 0.0, 0.0),
1595 ) -> None:
1596 """Rotate atoms via Euler angles (in degrees).
1598 See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation.
1600 Note that the rotations in this method are passive and applied **not**
1601 to the atomic coordinates in the present frame **but** the frame itself.
1603 Parameters
1604 ----------
1605 phi : float
1606 The 1st rotation angle around the z axis.
1607 theta : float
1608 Rotation around the x axis.
1609 psi : float
1610 2nd rotation around the z axis.
1611 center : Sequence[float], default = (0.0, 0.0, 0.0)
1612 The point to rotate about. A sequence of length 3 with the
1613 coordinates, or 'COM' to select the center of mass, 'COP' to
1614 select center of positions or 'COU' to select center of cell.
1616 """
1617 from scipy.spatial.transform import Rotation as R
1619 center = self._centering_as_array(center)
1621 # passive rotations (negative angles) for backward compatibility
1622 rotation = R.from_euler('zxz', (-phi, -theta, -psi), degrees=True)
1624 self.positions = rotation.apply(self.positions - center) + center
1626 def get_dihedral(self, a0, a1, a2, a3, mic=False):
1627 """Calculate dihedral angle.
1629 Calculate dihedral angle (in degrees) between the vectors a0->a1
1630 and a2->a3.
1632 Use mic=True to use the Minimum Image Convention and calculate the
1633 angle across periodic boundaries.
1634 """
1635 return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0]
1637 def get_dihedrals(self, indices, mic=False):
1638 """Calculate dihedral angles.
1640 Calculate dihedral angles (in degrees) between the list of vectors
1641 a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices.
1643 Use mic=True to use the Minimum Image Convention and calculate the
1644 angles across periodic boundaries.
1645 """
1646 from ase.geometry import get_dihedrals
1648 indices = np.array(indices)
1649 assert indices.shape[1] == 4
1651 a0s = self.positions[indices[:, 0]]
1652 a1s = self.positions[indices[:, 1]]
1653 a2s = self.positions[indices[:, 2]]
1654 a3s = self.positions[indices[:, 3]]
1656 # vectors 0->1, 1->2, 2->3
1657 v0 = a1s - a0s
1658 v1 = a2s - a1s
1659 v2 = a3s - a2s
1661 cell = None
1662 pbc = None
1664 if mic:
1665 cell = self.cell
1666 pbc = self.pbc
1668 return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc)
1670 def _masked_rotate(self, center, axis, diff, mask):
1671 # do rotation of subgroup by copying it to temporary atoms object
1672 # and then rotating that
1673 #
1674 # recursive object definition might not be the most elegant thing,
1675 # more generally useful might be a rotation function with a mask?
1676 group = self.__class__()
1677 for i in range(len(self)):
1678 if mask[i]:
1679 group += self[i]
1680 group.translate(-center)
1681 group.rotate(diff * 180 / pi, axis)
1682 group.translate(center)
1683 # set positions in original atoms object
1684 j = 0
1685 for i in range(len(self)):
1686 if mask[i]:
1687 self.positions[i] = group[j].position
1688 j += 1
1690 def set_dihedral(self, a1, a2, a3, a4, angle,
1691 mask=None, indices=None):
1692 """Set the dihedral angle (degrees) between vectors a1->a2 and
1693 a3->a4 by changing the atom indexed by a4.
1695 If mask is not None, all the atoms described in mask
1696 (read: the entire subgroup) are moved. Alternatively to the mask,
1697 the indices of the atoms to be rotated can be supplied. If both
1698 *mask* and *indices* are given, *indices* overwrites *mask*.
1700 **Important**: If *mask* or *indices* is given and does not contain
1701 *a4*, *a4* will NOT be moved. In most cases you therefore want
1702 to include *a4* in *mask*/*indices*.
1704 Example: the following defines a very crude
1705 ethane-like molecule and twists one half of it by 30 degrees.
1707 >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0],
1708 ... [1, 0, 0], [2, 1, 0], [2, -1, 0]])
1709 >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1])
1710 """
1712 angle *= pi / 180
1714 # if not provided, set mask to the last atom in the
1715 # dihedral description
1716 if mask is None and indices is None:
1717 mask = np.zeros(len(self))
1718 mask[a4] = 1
1719 elif indices is not None:
1720 mask = [index in indices for index in range(len(self))]
1722 # compute necessary in dihedral change, from current value
1723 current = self.get_dihedral(a1, a2, a3, a4) * pi / 180
1724 diff = angle - current
1725 axis = self.positions[a3] - self.positions[a2]
1726 center = self.positions[a3]
1727 self._masked_rotate(center, axis, diff, mask)
1729 def rotate_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None):
1730 """Rotate dihedral angle.
1732 Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a
1733 predefined dihedral angle, starting from its current configuration.
1734 """
1735 start = self.get_dihedral(a1, a2, a3, a4)
1736 self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices)
1738 def get_angle(self, a1, a2, a3, mic=False):
1739 """Get angle formed by three atoms.
1741 Calculate angle in degrees between the vectors a2->a1 and
1742 a2->a3.
1744 Use mic=True to use the Minimum Image Convention and calculate the
1745 angle across periodic boundaries.
1746 """
1747 return self.get_angles([[a1, a2, a3]], mic=mic)[0]
1749 def get_angles(self, indices, mic=False):
1750 """Get angle formed by three atoms for multiple groupings.
1752 Calculate angle in degrees between vectors between atoms a2->a1
1753 and a2->a3, where a1, a2, and a3 are in each row of indices.
1755 Use mic=True to use the Minimum Image Convention and calculate
1756 the angle across periodic boundaries.
1757 """
1758 from ase.geometry import get_angles
1760 indices = np.array(indices)
1761 assert indices.shape[1] == 3
1763 a1s = self.positions[indices[:, 0]]
1764 a2s = self.positions[indices[:, 1]]
1765 a3s = self.positions[indices[:, 2]]
1767 v12 = a1s - a2s
1768 v32 = a3s - a2s
1770 cell = None
1771 pbc = None
1773 if mic:
1774 cell = self.cell
1775 pbc = self.pbc
1777 return get_angles(v12, v32, cell=cell, pbc=pbc)
1779 def set_angle(self, a1, a2=None, a3=None, angle=None, mask=None,
1780 indices=None, add=False):
1781 """Set angle (in degrees) formed by three atoms.
1783 Sets the angle between vectors *a2*->*a1* and *a2*->*a3*.
1785 If *add* is `True`, the angle will be changed by the value given.
1787 Same usage as in :meth:`ase.Atoms.set_dihedral`.
1788 If *mask* and *indices*
1789 are given, *indices* overwrites *mask*. If *mask* and *indices*
1790 are not set, only *a3* is moved."""
1792 if any(a is None for a in [a2, a3, angle]):
1793 raise ValueError('a2, a3, and angle must not be None')
1795 # If not provided, set mask to the last atom in the angle description
1796 if mask is None and indices is None:
1797 mask = np.zeros(len(self))
1798 mask[a3] = 1
1799 elif indices is not None:
1800 mask = [index in indices for index in range(len(self))]
1802 if add:
1803 diff = angle
1804 else:
1805 # Compute necessary in angle change, from current value
1806 diff = angle - self.get_angle(a1, a2, a3)
1808 diff *= pi / 180
1809 # Do rotation of subgroup by copying it to temporary atoms object and
1810 # then rotating that
1811 v10 = self.positions[a1] - self.positions[a2]
1812 v12 = self.positions[a3] - self.positions[a2]
1813 v10 /= np.linalg.norm(v10)
1814 v12 /= np.linalg.norm(v12)
1815 axis = np.cross(v10, v12)
1816 center = self.positions[a2]
1817 self._masked_rotate(center, axis, diff, mask)
1819 def rattle(self, stdev=0.001, seed=None, rng=None):
1820 """Randomly displace atoms.
1822 This method adds random displacements to the atomic positions,
1823 taking a possible constraint into account. The random numbers are
1824 drawn from a normal distribution of standard deviation stdev.
1826 By default, the random number generator always uses the same seed (42)
1827 for repeatability. You can provide your own seed (an integer), or if you
1828 want the randomness to be different each time you run a script, then
1829 provide `rng=numpy.random`. For a parallel calculation, it is important
1830 to use the same seed on all processors! """
1832 if seed is not None and rng is not None:
1833 raise ValueError('Please do not provide both seed and rng.')
1835 if rng is None:
1836 if seed is None:
1837 seed = 42
1838 rng = np.random.RandomState(seed)
1839 positions = self.arrays['positions']
1840 self.set_positions(positions +
1841 rng.normal(scale=stdev, size=positions.shape))
1843 def get_distance(self, a0, a1, mic=False, vector=False):
1844 """Return distance between two atoms.
1846 Use mic=True to use the Minimum Image Convention.
1847 vector=True gives the distance vector (from a0 to a1).
1848 """
1849 return self.get_distances(a0, [a1], mic=mic, vector=vector)[0]
1851 def get_distances(self, a, indices, mic=False, vector=False):
1852 """Return distances of atom No.i with a list of atoms.
1854 Use mic=True to use the Minimum Image Convention.
1855 vector=True gives the distance vector (from a to self[indices]).
1856 """
1857 from ase.geometry import get_distances
1859 R = self.arrays['positions']
1860 p1 = [R[a]]
1861 p2 = R[indices]
1863 cell = None
1864 pbc = None
1866 if mic:
1867 cell = self.cell
1868 pbc = self.pbc
1870 D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc)
1872 if vector:
1873 D.shape = (-1, 3)
1874 return D
1875 else:
1876 D_len.shape = (-1,)
1877 return D_len
1879 def get_all_distances(self, mic=False, vector=False):
1880 """Return distances of all of the atoms with all of the atoms.
1882 Use mic=True to use the Minimum Image Convention.
1883 """
1884 from ase.geometry import get_distances
1886 R = self.arrays['positions']
1888 cell = None
1889 pbc = None
1891 if mic:
1892 cell = self.cell
1893 pbc = self.pbc
1895 D, D_len = get_distances(R, cell=cell, pbc=pbc)
1897 if vector:
1898 return D
1899 else:
1900 return D_len
1902 def set_distance(self, a0, a1, distance, fix=0.5, mic=False,
1903 mask=None, indices=None, add=False, factor=False):
1904 """Set the distance between two atoms.
1906 Set the distance between atoms *a0* and *a1* to *distance*.
1907 By default, the center of the two atoms will be fixed. Use
1908 *fix=0* to fix the first atom, *fix=1* to fix the second
1909 atom and *fix=0.5* (default) to fix the center of the bond.
1911 If *mask* or *indices* are set (*mask* overwrites *indices*),
1912 only the atoms defined there are moved
1913 (see :meth:`ase.Atoms.set_dihedral`).
1915 When *add* is true, the distance is changed by the value given.
1916 In combination
1917 with *factor* True, the value given is a factor scaling the distance.
1919 It is assumed that the atoms in *mask*/*indices* move together
1920 with *a1*. If *fix=1*, only *a0* will therefore be moved."""
1921 from ase.geometry import find_mic
1923 if a0 % len(self) == a1 % len(self):
1924 raise ValueError('a0 and a1 must not be the same')
1926 if add:
1927 oldDist = self.get_distance(a0, a1, mic=mic)
1928 if factor:
1929 newDist = oldDist * distance
1930 else:
1931 newDist = oldDist + distance
1932 self.set_distance(a0, a1, newDist, fix=fix, mic=mic,
1933 mask=mask, indices=indices, add=False,
1934 factor=False)
1935 return
1937 R = self.arrays['positions']
1938 D = np.array([R[a1] - R[a0]])
1940 if mic:
1941 D, D_len = find_mic(D, self.cell, self.pbc)
1942 else:
1943 D_len = np.array([np.sqrt((D**2).sum())])
1944 x = 1.0 - distance / D_len[0]
1946 if mask is None and indices is None:
1947 indices = [a0, a1]
1948 elif mask:
1949 indices = [i for i in range(len(self)) if mask[i]]
1951 for i in indices:
1952 if i == a0:
1953 R[a0] += (x * fix) * D[0]
1954 else:
1955 R[i] -= (x * (1.0 - fix)) * D[0]
1957 def get_scaled_positions(self, wrap=True):
1958 """Get positions relative to unit cell.
1960 If wrap is True, atoms outside the unit cell will be wrapped into
1961 the cell in those directions with periodic boundary conditions
1962 so that the scaled coordinates are between zero and one.
1964 If any cell vectors are zero, the corresponding coordinates
1965 are evaluated as if the cell were completed using
1966 ``cell.complete()``. This means coordinates will be Cartesian
1967 as long as the non-zero cell vectors span a Cartesian axis or
1968 plane."""
1970 fractional = self.cell.scaled_positions(self.positions)
1972 if wrap:
1973 for i, periodic in enumerate(self.pbc):
1974 if periodic:
1975 # Yes, we need to do it twice.
1976 # See the scaled_positions.py test.
1977 fractional[:, i] %= 1.0
1978 fractional[:, i] %= 1.0
1980 return fractional
1982 def set_scaled_positions(self, scaled):
1983 """Set positions relative to unit cell."""
1984 self.positions[:] = self.cell.cartesian_positions(scaled)
1986 def wrap(self, **wrap_kw):
1987 """Wrap positions to unit cell.
1989 Parameters:
1991 wrap_kw: (keyword=value) pairs
1992 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
1993 see :func:`ase.geometry.wrap_positions`
1994 """
1996 if 'pbc' not in wrap_kw:
1997 wrap_kw['pbc'] = self.pbc
1999 self.positions[:] = self.get_positions(wrap=True, **wrap_kw)
2001 def get_temperature(self):
2002 """Get the temperature in Kelvin."""
2003 ekin = self.get_kinetic_energy()
2004 return 2 * ekin / (self.get_number_of_degrees_of_freedom() * units.kB)
2006 def __eq__(self, other):
2007 """Check for identity of two atoms objects.
2009 Identity means: same positions, atomic numbers, unit cell and
2010 periodic boundary conditions."""
2011 if not isinstance(other, Atoms):
2012 return False
2013 a = self.arrays
2014 b = other.arrays
2015 return (len(self) == len(other) and
2016 (a['positions'] == b['positions']).all() and
2017 (a['numbers'] == b['numbers']).all() and
2018 (self.cell == other.cell).all() and
2019 (self.pbc == other.pbc).all())
2021 def __ne__(self, other):
2022 """Check if two atoms objects are not equal.
2024 Any differences in positions, atomic numbers, unit cell or
2025 periodic boundary condtions make atoms objects not equal.
2026 """
2027 eq = self.__eq__(other)
2028 if eq is NotImplemented:
2029 return eq
2030 else:
2031 return not eq
2033 # @deprecated('Please use atoms.cell.volume')
2034 # We kind of want to deprecate this, but the ValueError behaviour
2035 # might be desirable. Should we do this?
2036 def get_volume(self):
2037 """Get volume of unit cell."""
2038 if self.cell.rank != 3:
2039 raise ValueError(
2040 'You have {} lattice vectors: volume not defined'
2041 .format(self.cell.rank))
2042 return self.cell.volume
2044 @property
2045 def cell(self):
2046 """The :class:`ase.cell.Cell` for direct manipulation."""
2047 return self._cellobj
2049 @cell.setter
2050 def cell(self, cell):
2051 cell = Cell.ascell(cell)
2052 self._cellobj[:] = cell
2054 def write(self, filename, format=None, **kwargs):
2055 """Write atoms object to a file.
2057 see ase.io.write for formats.
2058 kwargs are passed to ase.io.write.
2059 """
2060 from ase.io import write
2061 write(filename, self, format, **kwargs)
2063 def iterimages(self):
2064 yield self
2066 def __ase_optimizable__(self):
2067 from ase.optimize.optimize import OptimizableAtoms
2068 return OptimizableAtoms(self)
2070 def edit(self):
2071 """Modify atoms interactively through ASE's GUI viewer.
2073 Conflicts leading to undesirable behaviour might arise
2074 when matplotlib has been pre-imported with certain
2075 incompatible backends and while trying to use the
2076 plot feature inside the interactive GUI. To circumvent,
2077 please set matplotlib.use('gtk') before calling this
2078 method.
2079 """
2080 from ase.gui.gui import GUI
2081 from ase.gui.images import Images
2082 images = Images([self])
2083 gui = GUI(images)
2084 gui.run()
2087def string2vector(v):
2088 if isinstance(v, str):
2089 if v[0] == '-':
2090 return -string2vector(v[1:])
2091 w = np.zeros(3)
2092 w['xyz'.index(v)] = 1.0
2093 return w
2094 return np.array(v, float)
2097def default(data, dflt):
2098 """Helper function for setting default values."""
2099 if data is None:
2100 return None
2101 elif isinstance(data, (list, tuple)):
2102 newdata = []
2103 allnone = True
2104 for x in data:
2105 if x is None:
2106 newdata.append(dflt)
2107 else:
2108 newdata.append(x)
2109 allnone = False
2110 if allnone:
2111 return None
2112 return newdata
2113 else:
2114 return data