Coverage for /builds/debichem-team/python-ase/ase/atoms.py: 93.66%
978 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1# Copyright 2008, 2009 CAMd
2# (see accompanying license files for details).
4"""Definition of the Atoms class.
6This module defines the central object in the ASE package: the Atoms
7object.
8"""
9import copy
10import numbers
11from math import cos, pi, sin
13import numpy as np
15import ase.units as units
16from ase.atom import Atom
17from ase.cell import Cell
18from ase.data import atomic_masses, atomic_masses_common
19from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
20from ase.symbols import Symbols, symbols2numbers
21from ase.utils import deprecated, string2index
24class Atoms:
25 """Atoms object.
27 The Atoms object can represent an isolated molecule, or a
28 periodically repeated structure. It has a unit cell and
29 there may be periodic boundary conditions along any of the three
30 unit cell axes.
31 Information about the atoms (atomic numbers and position) is
32 stored in ndarrays. Optionally, there can be information about
33 tags, momenta, masses, magnetic moments and charges.
35 In order to calculate energies, forces and stresses, a calculator
36 object has to attached to the atoms object.
38 Parameters
39 ----------
40 symbols : str | list[str] | list[Atom]
41 Chemical formula, a list of chemical symbols, or list of
42 :class:`~ase.Atom` objects (mutually exclusive with ``numbers``).
44 - ``'H2O'``
45 - ``'COPt12'``
46 - ``['H', 'H', 'O']``
47 - ``[Atom('Ne', (x, y, z)), ...]``
49 positions : list[tuple[float, float, float]]
50 Atomic positions in Cartesian coordinates
51 (mutually exclusive with ``scaled_positions``).
52 Anything that can be converted to an ndarray of shape (n, 3) works:
53 [(x0, y0, z0), (x1, y1, z1), ...].
54 scaled_positions : list[tuple[float, float, float]]
55 Atomic positions in units of the unit cell
56 (mutually exclusive with ``positions``).
57 numbers : list[int]
58 Atomic numbers (mutually exclusive with ``symbols``).
59 tags : list[int]
60 Special purpose tags.
61 momenta : list[tuple[float, float, float]]
62 Momenta for all atoms in Cartesian coordinates
63 (mutually exclusive with ``velocities``).
64 velocities : list[tuple[float, float, float]]
65 Velocities for all atoms in Cartesian coordinates
66 (mutually exclusive with ``momenta``).
67 masses : list[float]
68 Atomic masses in atomic units.
69 magmoms : list[float] | list[tuple[float, float, float]]
70 Magnetic moments. Can be either a single value for each atom
71 for collinear calculations or three numbers for each atom for
72 non-collinear calculations.
73 charges : list[float]
74 Initial atomic charges.
75 cell : 3x3 matrix or length 3 or 6 vector, default: (0, 0, 0)
76 Unit cell vectors. Can also be given as just three
77 numbers for orthorhombic cells, or 6 numbers, where
78 first three are lengths of unit cell vectors, and the
79 other three are angles between them (in degrees), in following order:
80 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)].
81 First vector will lie in x-direction, second in xy-plane,
82 and the third one in z-positive subspace.
83 celldisp : tuple[float, float, float], default: (0, 0, 0)
84 Unit cell displacement vector. To visualize a displaced cell
85 around the center of mass of a Systems of atoms.
86 pbc : bool | tuple[bool, bool, bool], default: False
87 Periodic boundary conditions flags.
89 - ``True``
90 - ``False``
91 - ``0``
92 - ``1``
93 - ``(1, 1, 0)``
94 - ``(True, False, False)``
96 constraint : constraint object(s)
97 One or more ASE constraints applied during structure optimization.
98 calculator : calculator object
99 ASE calculator to obtain energies and atomic forces.
100 info : dict | None, default: None
101 Dictionary with additional information about the system.
102 The following keys may be used by ASE:
104 - spacegroup: :class:`~ase.spacegroup.Spacegroup` instance
105 - unit_cell: 'conventional' | 'primitive' | int | 3 ints
106 - adsorbate_info: Information about special adsorption sites
108 Items in the info attribute survives copy and slicing and can
109 be stored in and retrieved from trajectory files given that the
110 key is a string, the value is JSON-compatible and, if the value is a
111 user-defined object, its base class is importable. One should
112 not make any assumptions about the existence of keys.
114 Examples
115 --------
116 >>> from ase import Atom
118 N2 molecule (These three are equivalent):
120 >>> d = 1.104 # N2 bondlength
121 >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)])
122 >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)])
123 >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))])
125 FCC gold:
127 >>> a = 4.05 # Gold lattice constant
128 >>> b = a / 2
129 >>> fcc = Atoms('Au',
130 ... cell=[(0, b, b), (b, 0, b), (b, b, 0)],
131 ... pbc=True)
133 Hydrogen wire:
135 >>> d = 0.9 # H-H distance
136 >>> h = Atoms('H', positions=[(0, 0, 0)],
137 ... cell=(d, 0, 0),
138 ... pbc=(1, 0, 0))
139 """
141 ase_objtype = 'atoms' # For JSONability
143 def __init__(self, symbols=None,
144 positions=None, numbers=None,
145 tags=None, momenta=None, masses=None,
146 magmoms=None, charges=None,
147 scaled_positions=None,
148 cell=None, pbc=None, celldisp=None,
149 constraint=None,
150 calculator=None,
151 info=None,
152 velocities=None):
154 self._cellobj = Cell.new()
155 self._pbc = np.zeros(3, bool)
157 atoms = None
159 if hasattr(symbols, 'get_positions'):
160 atoms = symbols
161 symbols = None
162 elif (isinstance(symbols, (list, tuple)) and
163 len(symbols) > 0 and isinstance(symbols[0], Atom)):
164 # Get data from a list or tuple of Atom objects:
165 data = [[atom.get_raw(name) for atom in symbols]
166 for name in
167 ['position', 'number', 'tag', 'momentum',
168 'mass', 'magmom', 'charge']]
169 atoms = self.__class__(None, *data)
170 symbols = None
172 if atoms is not None:
173 # Get data from another Atoms object:
174 if scaled_positions is not None:
175 raise NotImplementedError
176 if symbols is None and numbers is None:
177 numbers = atoms.get_atomic_numbers()
178 if positions is None:
179 positions = atoms.get_positions()
180 if tags is None and atoms.has('tags'):
181 tags = atoms.get_tags()
182 if momenta is None and atoms.has('momenta'):
183 momenta = atoms.get_momenta()
184 if magmoms is None and atoms.has('initial_magmoms'):
185 magmoms = atoms.get_initial_magnetic_moments()
186 if masses is None and atoms.has('masses'):
187 masses = atoms.get_masses()
188 if charges is None and atoms.has('initial_charges'):
189 charges = atoms.get_initial_charges()
190 if cell is None:
191 cell = atoms.get_cell()
192 if celldisp is None:
193 celldisp = atoms.get_celldisp()
194 if pbc is None:
195 pbc = atoms.get_pbc()
196 if constraint is None:
197 constraint = [c.copy() for c in atoms.constraints]
198 if calculator is None:
199 calculator = atoms.calc
200 if info is None:
201 info = copy.deepcopy(atoms.info)
203 self.arrays = {}
205 if symbols is not None and numbers is not None:
206 raise TypeError('Use only one of "symbols" and "numbers".')
207 if symbols is not None:
208 numbers = symbols2numbers(symbols)
209 elif numbers is None:
210 if positions is not None:
211 natoms = len(positions)
212 elif scaled_positions is not None:
213 natoms = len(scaled_positions)
214 else:
215 natoms = 0
216 numbers = np.zeros(natoms, int)
217 self.new_array('numbers', numbers, int)
219 if self.numbers.ndim != 1:
220 raise ValueError('"numbers" must be 1-dimensional.')
222 if cell is None:
223 cell = np.zeros((3, 3))
224 self.set_cell(cell)
226 if celldisp is None:
227 celldisp = np.zeros(shape=(3, 1))
228 self.set_celldisp(celldisp)
230 if positions is None:
231 if scaled_positions is None:
232 positions = np.zeros((len(self.arrays['numbers']), 3))
233 else:
234 assert self.cell.rank == 3
235 positions = np.dot(scaled_positions, self.cell)
236 else:
237 if scaled_positions is not None:
238 raise TypeError(
239 'Use only one of "positions" and "scaled_positions".')
240 self.new_array('positions', positions, float, (3,))
242 self.set_constraint(constraint)
243 self.set_tags(default(tags, 0))
244 self.set_masses(default(masses, None))
245 self.set_initial_magnetic_moments(default(magmoms, 0.0))
246 self.set_initial_charges(default(charges, 0.0))
247 if pbc is None:
248 pbc = False
249 self.set_pbc(pbc)
250 self.set_momenta(default(momenta, (0.0, 0.0, 0.0)),
251 apply_constraint=False)
253 if velocities is not None:
254 if momenta is None:
255 self.set_velocities(velocities)
256 else:
257 raise TypeError(
258 'Use only one of "momenta" and "velocities".')
260 if info is None:
261 self.info = {}
262 else:
263 self.info = dict(info)
265 self.calc = calculator
267 @property
268 def symbols(self):
269 """Get chemical symbols as a :class:`ase.symbols.Symbols` object.
271 The object works like ``atoms.numbers`` except its values
272 are strings. It supports in-place editing."""
273 return Symbols(self.numbers)
275 @symbols.setter
276 def symbols(self, obj):
277 new_symbols = Symbols.fromsymbols(obj)
278 self.numbers[:] = new_symbols.numbers
280 @deprecated("Please use atoms.calc = calc", FutureWarning)
281 def set_calculator(self, calc=None):
282 """Attach calculator object.
284 .. deprecated:: 3.20.0
285 Please use the equivalent ``atoms.calc = calc`` instead of this
286 method.
287 """
289 self.calc = calc
291 @deprecated("Please use atoms.calc", FutureWarning)
292 def get_calculator(self):
293 """Get currently attached calculator object.
295 .. deprecated:: 3.20.0
296 Please use the equivalent ``atoms.calc`` instead of
297 ``atoms.get_calculator()``.
298 """
300 return self.calc
302 @property
303 def calc(self):
304 """Calculator object."""
305 return self._calc
307 @calc.setter
308 def calc(self, calc):
309 self._calc = calc
310 if hasattr(calc, 'set_atoms'):
311 calc.set_atoms(self)
313 @calc.deleter
314 @deprecated('Please use atoms.calc = None', FutureWarning)
315 def calc(self):
316 """Delete calculator
318 .. deprecated:: 3.20.0
319 Please use ``atoms.calc = None``
320 """
321 self._calc = None
323 @property
324 @deprecated('Please use atoms.cell.rank instead', DeprecationWarning)
325 def number_of_lattice_vectors(self):
326 """Number of (non-zero) lattice vectors.
328 .. deprecated:: 3.21.0
329 Please use ``atoms.cell.rank`` instead
330 """
331 return self.cell.rank
333 def set_constraint(self, constraint=None):
334 """Apply one or more constrains.
336 The *constraint* argument must be one constraint object or a
337 list of constraint objects."""
338 if constraint is None:
339 self._constraints = []
340 else:
341 if isinstance(constraint, list):
342 self._constraints = constraint
343 elif isinstance(constraint, tuple):
344 self._constraints = list(constraint)
345 else:
346 self._constraints = [constraint]
348 def _get_constraints(self):
349 return self._constraints
351 def _del_constraints(self):
352 self._constraints = []
354 constraints = property(_get_constraints, set_constraint, _del_constraints,
355 'Constraints of the atoms.')
357 def get_number_of_degrees_of_freedom(self):
358 """Calculate the number of degrees of freedom in the system."""
359 return len(self) * 3 - sum(
360 c.get_removed_dof(self) for c in self._constraints
361 )
363 def set_cell(self, cell, scale_atoms=False, apply_constraint=True):
364 """Set unit cell vectors.
366 Parameters:
368 cell: 3x3 matrix or length 3 or 6 vector
369 Unit cell. A 3x3 matrix (the three unit cell vectors) or
370 just three numbers for an orthorhombic cell. Another option is
371 6 numbers, which describes unit cell with lengths of unit cell
372 vectors and with angles between them (in degrees), in following
373 order: [len(a), len(b), len(c), angle(b,c), angle(a,c),
374 angle(a,b)]. First vector will lie in x-direction, second in
375 xy-plane, and the third one in z-positive subspace.
376 scale_atoms: bool
377 Fix atomic positions or move atoms with the unit cell?
378 Default behavior is to *not* move the atoms (scale_atoms=False).
379 apply_constraint: bool
380 Whether to apply constraints to the given cell.
382 Examples:
384 Two equivalent ways to define an orthorhombic cell:
386 >>> atoms = Atoms('He')
387 >>> a, b, c = 7, 7.5, 8
388 >>> atoms.set_cell([a, b, c])
389 >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)])
391 FCC unit cell:
393 >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)])
395 Hexagonal unit cell:
397 >>> atoms.set_cell([a, a, c, 90, 90, 120])
399 Rhombohedral unit cell:
401 >>> alpha = 77
402 >>> atoms.set_cell([a, a, a, alpha, alpha, alpha])
403 """
405 # Override pbcs if and only if given a Cell object:
406 cell = Cell.new(cell)
408 # XXX not working well during initialize due to missing _constraints
409 if apply_constraint and hasattr(self, '_constraints'):
410 for constraint in self.constraints:
411 if hasattr(constraint, 'adjust_cell'):
412 constraint.adjust_cell(self, cell)
414 if scale_atoms:
415 M = np.linalg.solve(self.cell.complete(), cell.complete())
416 self.positions[:] = np.dot(self.positions, M)
418 self.cell[:] = cell
420 def set_celldisp(self, celldisp):
421 """Set the unit cell displacement vectors."""
422 celldisp = np.array(celldisp, float)
423 self._celldisp = celldisp
425 def get_celldisp(self):
426 """Get the unit cell displacement vectors."""
427 return self._celldisp.copy()
429 def get_cell(self, complete=False):
430 """Get the three unit cell vectors as a `class`:ase.cell.Cell` object.
432 The Cell object resembles a 3x3 ndarray, and cell[i, j]
433 is the jth Cartesian coordinate of the ith cell vector."""
434 if complete:
435 cell = self.cell.complete()
436 else:
437 cell = self.cell.copy()
439 return cell
441 @deprecated('Please use atoms.cell.cellpar() instead', DeprecationWarning)
442 def get_cell_lengths_and_angles(self):
443 """Get unit cell parameters. Sequence of 6 numbers.
445 First three are unit cell vector lengths and second three
446 are angles between them::
448 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
450 in degrees.
452 .. deprecated:: 3.21.0
453 Please use ``atoms.cell.cellpar()`` instead
454 """
455 return self.cell.cellpar()
457 @deprecated('Please use atoms.cell.reciprocal()', DeprecationWarning)
458 def get_reciprocal_cell(self):
459 """Get the three reciprocal lattice vectors as a 3x3 ndarray.
461 Note that the commonly used factor of 2 pi for Fourier
462 transforms is not included here.
464 .. deprecated:: 3.21.0
465 Please use ``atoms.cell.reciprocal()``
466 """
467 return self.cell.reciprocal()
469 @property
470 def pbc(self):
471 """Reference to pbc-flags for in-place manipulations."""
472 return self._pbc
474 @pbc.setter
475 def pbc(self, pbc):
476 self._pbc[:] = pbc
478 def set_pbc(self, pbc):
479 """Set periodic boundary condition flags."""
480 self.pbc = pbc
482 def get_pbc(self):
483 """Get periodic boundary condition flags."""
484 return self.pbc.copy()
486 def new_array(self, name, a, dtype=None, shape=None):
487 """Add new array.
489 If *shape* is not *None*, the shape of *a* will be checked."""
491 if dtype is not None:
492 a = np.array(a, dtype, order='C')
493 if len(a) == 0 and shape is not None:
494 a.shape = (-1,) + shape
495 else:
496 if not a.flags['C_CONTIGUOUS']:
497 a = np.ascontiguousarray(a)
498 else:
499 a = a.copy()
501 if name in self.arrays:
502 raise RuntimeError(f'Array {name} already present')
504 for b in self.arrays.values():
505 if len(a) != len(b):
506 raise ValueError('Array "%s" has wrong length: %d != %d.' %
507 (name, len(a), len(b)))
508 break
510 if shape is not None and a.shape[1:] != shape:
511 raise ValueError(
512 f'Array "{name}" has wrong shape {a.shape} != '
513 f'{(a.shape[0:1] + shape)}.')
515 self.arrays[name] = a
517 def get_array(self, name, copy=True):
518 """Get an array.
520 Returns a copy unless the optional argument copy is false.
521 """
522 if copy:
523 return self.arrays[name].copy()
524 else:
525 return self.arrays[name]
527 def set_array(self, name, a, dtype=None, shape=None):
528 """Update array.
530 If *shape* is not *None*, the shape of *a* will be checked.
531 If *a* is *None*, then the array is deleted."""
533 b = self.arrays.get(name)
534 if b is None:
535 if a is not None:
536 self.new_array(name, a, dtype, shape)
537 else:
538 if a is None:
539 del self.arrays[name]
540 else:
541 a = np.asarray(a)
542 if a.shape != b.shape:
543 raise ValueError(
544 f'Array "{name}" has wrong shape '
545 f'{a.shape} != {b.shape}.')
546 b[:] = a
548 def has(self, name):
549 """Check for existence of array.
551 name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms',
552 'initial_charges'."""
553 # XXX extend has to calculator properties
554 return name in self.arrays
556 def set_atomic_numbers(self, numbers):
557 """Set atomic numbers."""
558 self.set_array('numbers', numbers, int, ())
560 def get_atomic_numbers(self):
561 """Get integer array of atomic numbers."""
562 return self.arrays['numbers'].copy()
564 def get_chemical_symbols(self):
565 """Get list of chemical symbol strings.
567 Equivalent to ``list(atoms.symbols)``."""
568 return list(self.symbols)
570 def set_chemical_symbols(self, symbols):
571 """Set chemical symbols."""
572 self.set_array('numbers', symbols2numbers(symbols), int, ())
574 def get_chemical_formula(self, mode='hill', empirical=False):
575 """Get the chemical formula as a string based on the chemical symbols.
577 Parameters:
579 mode: str
580 There are four different modes available:
582 'all': The list of chemical symbols are contracted to a string,
583 e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'.
585 'reduce': The same as 'all' where repeated elements are contracted
586 to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to
587 'CH3OCH3'.
589 'hill': The list of chemical symbols are contracted to a string
590 following the Hill notation (alphabetical order with C and H
591 first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to
592 'H2O4S'. This is default.
594 'metal': The list of chemical symbols (alphabetical metals,
595 and alphabetical non-metals)
597 empirical, bool (optional, default=False)
598 Divide the symbol counts by their greatest common divisor to yield
599 an empirical formula. Only for mode `metal` and `hill`.
600 """
601 return self.symbols.get_chemical_formula(mode, empirical)
603 def set_tags(self, tags):
604 """Set tags for all atoms. If only one tag is supplied, it is
605 applied to all atoms."""
606 if isinstance(tags, int):
607 tags = [tags] * len(self)
608 self.set_array('tags', tags, int, ())
610 def get_tags(self):
611 """Get integer array of tags."""
612 if 'tags' in self.arrays:
613 return self.arrays['tags'].copy()
614 else:
615 return np.zeros(len(self), int)
617 def set_momenta(self, momenta, apply_constraint=True):
618 """Set momenta."""
619 if (apply_constraint and len(self.constraints) > 0 and
620 momenta is not None):
621 momenta = np.array(momenta) # modify a copy
622 for constraint in self.constraints:
623 if hasattr(constraint, 'adjust_momenta'):
624 constraint.adjust_momenta(self, momenta)
625 self.set_array('momenta', momenta, float, (3,))
627 def set_velocities(self, velocities):
628 """Set the momenta by specifying the velocities."""
629 self.set_momenta(self.get_masses()[:, np.newaxis] * velocities)
631 def get_momenta(self):
632 """Get array of momenta."""
633 if 'momenta' in self.arrays:
634 return self.arrays['momenta'].copy()
635 else:
636 return np.zeros((len(self), 3))
638 def set_masses(self, masses='defaults'):
639 """Set atomic masses in atomic mass units.
641 The array masses should contain a list of masses. In case
642 the masses argument is not given or for those elements of the
643 masses list that are None, standard values are set."""
645 if isinstance(masses, str):
646 if masses == 'defaults':
647 masses = atomic_masses[self.arrays['numbers']]
648 elif masses == 'most_common':
649 masses = atomic_masses_common[self.arrays['numbers']]
650 elif masses is None:
651 pass
652 elif not isinstance(masses, np.ndarray):
653 masses = list(masses)
654 for i, mass in enumerate(masses):
655 if mass is None:
656 masses[i] = atomic_masses[self.numbers[i]]
657 self.set_array('masses', masses, float, ())
659 def get_masses(self):
660 """Get array of masses in atomic mass units."""
661 if 'masses' in self.arrays:
662 return self.arrays['masses'].copy()
663 else:
664 return atomic_masses[self.arrays['numbers']]
666 def set_initial_magnetic_moments(self, magmoms=None):
667 """Set the initial magnetic moments.
669 Use either one or three numbers for every atom (collinear
670 or non-collinear spins)."""
672 if magmoms is None:
673 self.set_array('initial_magmoms', None)
674 else:
675 magmoms = np.asarray(magmoms)
676 self.set_array('initial_magmoms', magmoms, float,
677 magmoms.shape[1:])
679 def get_initial_magnetic_moments(self):
680 """Get array of initial magnetic moments."""
681 if 'initial_magmoms' in self.arrays:
682 return self.arrays['initial_magmoms'].copy()
683 else:
684 return np.zeros(len(self))
686 def get_magnetic_moments(self):
687 """Get calculated local magnetic moments."""
688 if self._calc is None:
689 raise RuntimeError('Atoms object has no calculator.')
690 return self._calc.get_magnetic_moments(self)
692 def get_magnetic_moment(self):
693 """Get calculated total magnetic moment."""
694 if self._calc is None:
695 raise RuntimeError('Atoms object has no calculator.')
696 return self._calc.get_magnetic_moment(self)
698 def set_initial_charges(self, charges=None):
699 """Set the initial charges."""
701 if charges is None:
702 self.set_array('initial_charges', None)
703 else:
704 self.set_array('initial_charges', charges, float, ())
706 def get_initial_charges(self):
707 """Get array of initial charges."""
708 if 'initial_charges' in self.arrays:
709 return self.arrays['initial_charges'].copy()
710 else:
711 return np.zeros(len(self))
713 def get_charges(self):
714 """Get calculated charges."""
715 if self._calc is None:
716 raise RuntimeError('Atoms object has no calculator.')
717 try:
718 return self._calc.get_charges(self)
719 except AttributeError:
720 from ase.calculators.calculator import PropertyNotImplementedError
721 raise PropertyNotImplementedError
723 def set_positions(self, newpositions, apply_constraint=True):
724 """Set positions, honoring any constraints. To ignore constraints,
725 use *apply_constraint=False*."""
726 if self.constraints and apply_constraint:
727 newpositions = np.array(newpositions, float)
728 for constraint in self.constraints:
729 constraint.adjust_positions(self, newpositions)
731 self.set_array('positions', newpositions, shape=(3,))
733 def get_positions(self, wrap=False, **wrap_kw):
734 """Get array of positions.
736 Parameters:
738 wrap: bool
739 wrap atoms back to the cell before returning positions
740 wrap_kw: (keyword=value) pairs
741 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
742 see :func:`ase.geometry.wrap_positions`
743 """
744 from ase.geometry import wrap_positions
745 if wrap:
746 if 'pbc' not in wrap_kw:
747 wrap_kw['pbc'] = self.pbc
748 return wrap_positions(self.positions, self.cell, **wrap_kw)
749 else:
750 return self.arrays['positions'].copy()
752 def get_potential_energy(self, force_consistent=False,
753 apply_constraint=True):
754 """Calculate potential energy.
756 Ask the attached calculator to calculate the potential energy and
757 apply constraints. Use *apply_constraint=False* to get the raw
758 forces.
760 When supported by the calculator, either the energy extrapolated
761 to zero Kelvin or the energy consistent with the forces (the free
762 energy) can be returned.
763 """
764 if self._calc is None:
765 raise RuntimeError('Atoms object has no calculator.')
766 if force_consistent:
767 energy = self._calc.get_potential_energy(
768 self, force_consistent=force_consistent)
769 else:
770 energy = self._calc.get_potential_energy(self)
771 if apply_constraint:
772 for constraint in self.constraints:
773 if hasattr(constraint, 'adjust_potential_energy'):
774 energy += constraint.adjust_potential_energy(self)
775 return energy
777 def get_properties(self, properties):
778 """This method is experimental; currently for internal use."""
779 # XXX Something about constraints.
780 if self._calc is None:
781 raise RuntimeError('Atoms object has no calculator.')
782 return self._calc.calculate_properties(self, properties)
784 def get_potential_energies(self):
785 """Calculate the potential energies of all the atoms.
787 Only available with calculators supporting per-atom energies
788 (e.g. classical potentials).
789 """
790 if self._calc is None:
791 raise RuntimeError('Atoms object has no calculator.')
792 return self._calc.get_potential_energies(self)
794 def get_kinetic_energy(self):
795 """Get the kinetic energy."""
796 momenta = self.arrays.get('momenta')
797 if momenta is None:
798 return 0.0
799 return 0.5 * np.vdot(momenta, self.get_velocities())
801 def get_velocities(self):
802 """Get array of velocities."""
803 momenta = self.get_momenta()
804 masses = self.get_masses()
805 return momenta / masses[:, np.newaxis]
807 def get_total_energy(self):
808 """Get the total energy - potential plus kinetic energy."""
809 return self.get_potential_energy() + self.get_kinetic_energy()
811 def get_forces(self, apply_constraint=True, md=False):
812 """Calculate atomic forces.
814 Ask the attached calculator to calculate the forces and apply
815 constraints. Use *apply_constraint=False* to get the raw
816 forces.
818 For molecular dynamics (md=True) we don't apply the constraint
819 to the forces but to the momenta. When holonomic constraints for
820 rigid linear triatomic molecules are present, ask the constraints
821 to redistribute the forces within each triple defined in the
822 constraints (required for molecular dynamics with this type of
823 constraints)."""
825 if self._calc is None:
826 raise RuntimeError('Atoms object has no calculator.')
827 forces = self._calc.get_forces(self)
829 if apply_constraint:
830 # We need a special md flag here because for MD we want
831 # to skip real constraints but include special "constraints"
832 # Like Hookean.
833 for constraint in self.constraints:
834 if md and hasattr(constraint, 'redistribute_forces_md'):
835 constraint.redistribute_forces_md(self, forces)
836 if not md or hasattr(constraint, 'adjust_potential_energy'):
837 constraint.adjust_forces(self, forces)
838 return forces
840 # Informs calculators (e.g. Asap) that ideal gas contribution is added here.
841 _ase_handles_dynamic_stress = True
843 def get_stress(self, voigt=True, apply_constraint=True,
844 include_ideal_gas=False):
845 """Calculate stress tensor.
847 Returns an array of the six independent components of the
848 symmetric stress tensor, in the traditional Voigt order
849 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt
850 order.
852 The ideal gas contribution to the stresses is added if the
853 atoms have momenta and ``include_ideal_gas`` is set to True.
854 """
856 if self._calc is None:
857 raise RuntimeError('Atoms object has no calculator.')
859 stress = self._calc.get_stress(self)
860 shape = stress.shape
862 if shape == (3, 3):
863 # Convert to the Voigt form before possibly applying
864 # constraints and adding the dynamic part of the stress
865 # (the "ideal gas contribution").
866 stress = full_3x3_to_voigt_6_stress(stress)
867 else:
868 assert shape == (6,)
870 if apply_constraint:
871 for constraint in self.constraints:
872 if hasattr(constraint, 'adjust_stress'):
873 constraint.adjust_stress(self, stress)
875 # Add ideal gas contribution, if applicable
876 if include_ideal_gas and self.has('momenta'):
877 stress += self.get_kinetic_stress()
879 if voigt:
880 return stress
881 else:
882 return voigt_6_to_full_3x3_stress(stress)
884 def get_stresses(self, include_ideal_gas=False, voigt=True):
885 """Calculate the stress-tensor of all the atoms.
887 Only available with calculators supporting per-atom energies and
888 stresses (e.g. classical potentials). Even for such calculators
889 there is a certain arbitrariness in defining per-atom stresses.
891 The ideal gas contribution to the stresses is added if the
892 atoms have momenta and ``include_ideal_gas`` is set to True.
893 """
894 if self._calc is None:
895 raise RuntimeError('Atoms object has no calculator.')
896 stresses = self._calc.get_stresses(self)
898 # make sure `stresses` are in voigt form
899 if np.shape(stresses)[1:] == (3, 3):
900 stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses]
901 stresses = np.array(stresses_voigt)
903 # REMARK: The ideal gas contribution is intensive, i.e., the volume
904 # is divided out. We currently don't check if `stresses` are intensive
905 # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`.
906 # It might be good to check this here, but adds computational overhead.
908 if include_ideal_gas and self.has('momenta'):
909 stresses += self.get_kinetic_stresses()
911 if voigt:
912 return stresses
913 else:
914 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
915 return np.array(stresses_3x3)
917 def get_kinetic_stress(self, voigt=True):
918 """Calculate the kinetic part of the Virial stress tensor."""
919 stress = np.zeros(6) # Voigt notation
920 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
921 p = self.get_momenta()
922 masses = self.get_masses()
923 invmass = 1.0 / masses
924 invvol = 1.0 / self.get_volume()
925 for alpha in range(3):
926 for beta in range(alpha, 3):
927 stress[stresscomp[alpha, beta]] -= (
928 p[:, alpha] * p[:, beta] * invmass).sum() * invvol
930 if voigt:
931 return stress
932 else:
933 return voigt_6_to_full_3x3_stress(stress)
935 def get_kinetic_stresses(self, voigt=True):
936 """Calculate the kinetic part of the Virial stress of all the atoms."""
937 stresses = np.zeros((len(self), 6)) # Voigt notation
938 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
939 if hasattr(self._calc, 'get_atomic_volumes'):
940 invvol = 1.0 / self._calc.get_atomic_volumes()
941 else:
942 invvol = self.get_global_number_of_atoms() / self.get_volume()
943 p = self.get_momenta()
944 invmass = 1.0 / self.get_masses()
945 for alpha in range(3):
946 for beta in range(alpha, 3):
947 stresses[:, stresscomp[alpha, beta]] -= (
948 p[:, alpha] * p[:, beta] * invmass * invvol)
950 if voigt:
951 return stresses
952 else:
953 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
954 return np.array(stresses_3x3)
956 def get_dipole_moment(self):
957 """Calculate the electric dipole moment for the atoms object.
959 Only available for calculators which has a get_dipole_moment()
960 method."""
962 if self._calc is None:
963 raise RuntimeError('Atoms object has no calculator.')
964 return self._calc.get_dipole_moment(self)
966 def copy(self):
967 """Return a copy."""
968 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
969 celldisp=self._celldisp.copy())
971 atoms.arrays = {}
972 for name, a in self.arrays.items():
973 atoms.arrays[name] = a.copy()
974 atoms.constraints = copy.deepcopy(self.constraints)
975 return atoms
977 def todict(self):
978 """For basic JSON (non-database) support."""
979 d = dict(self.arrays)
980 d['cell'] = np.asarray(self.cell)
981 d['pbc'] = self.pbc
982 if self._celldisp.any():
983 d['celldisp'] = self._celldisp
984 if self.constraints:
985 d['constraints'] = self.constraints
986 if self.info:
987 d['info'] = self.info
988 # Calculator... trouble.
989 return d
991 @classmethod
992 def fromdict(cls, dct):
993 """Rebuild atoms object from dictionary representation (todict)."""
994 dct = dct.copy()
995 kw = {name: dct.pop(name)
996 for name in ['numbers', 'positions', 'cell', 'pbc']}
997 constraints = dct.pop('constraints', None)
998 if constraints:
999 from ase.constraints import dict2constraint
1000 constraints = [dict2constraint(d) for d in constraints]
1002 info = dct.pop('info', None)
1004 atoms = cls(constraint=constraints,
1005 celldisp=dct.pop('celldisp', None),
1006 info=info, **kw)
1007 natoms = len(atoms)
1009 # Some arrays are named differently from the atoms __init__ keywords.
1010 # Also, there may be custom arrays. Hence we set them directly:
1011 for name, arr in dct.items():
1012 assert len(arr) == natoms, name
1013 assert isinstance(arr, np.ndarray)
1014 atoms.arrays[name] = arr
1015 return atoms
1017 def __len__(self):
1018 return len(self.arrays['positions'])
1020 @deprecated(
1021 "Please use len(self) or, if your atoms are distributed, "
1022 "self.get_global_number_of_atoms.",
1023 category=FutureWarning,
1024 )
1025 def get_number_of_atoms(self):
1026 """
1027 .. deprecated:: 3.18.1
1028 You probably want ``len(atoms)``. Or if your atoms are distributed,
1029 use (and see) :func:`get_global_number_of_atoms()`.
1030 """
1031 return len(self)
1033 def get_global_number_of_atoms(self):
1034 """Returns the global number of atoms in a distributed-atoms parallel
1035 simulation.
1037 DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING!
1039 Equivalent to len(atoms) in the standard ASE Atoms class. You should
1040 normally use len(atoms) instead. This function's only purpose is to
1041 make compatibility between ASE and Asap easier to maintain by having a
1042 few places in ASE use this function instead. It is typically only
1043 when counting the global number of degrees of freedom or in similar
1044 situations.
1045 """
1046 return len(self)
1048 def __repr__(self):
1049 tokens = []
1051 N = len(self)
1052 if N <= 60:
1053 symbols = self.get_chemical_formula('reduce')
1054 else:
1055 symbols = self.get_chemical_formula('hill')
1056 tokens.append(f"symbols='{symbols}'")
1058 if self.pbc.any() and not self.pbc.all():
1059 tokens.append(f'pbc={self.pbc.tolist()}')
1060 else:
1061 tokens.append(f'pbc={self.pbc[0]}')
1063 cell = self.cell
1064 if cell:
1065 if cell.orthorhombic:
1066 cell = cell.lengths().tolist()
1067 else:
1068 cell = cell.tolist()
1069 tokens.append(f'cell={cell}')
1071 for name in sorted(self.arrays):
1072 if name in ['numbers', 'positions']:
1073 continue
1074 tokens.append(f'{name}=...')
1076 if self.constraints:
1077 if len(self.constraints) == 1:
1078 constraint = self.constraints[0]
1079 else:
1080 constraint = self.constraints
1081 tokens.append(f'constraint={constraint!r}')
1083 if self._calc is not None:
1084 tokens.append('calculator={}(...)'
1085 .format(self._calc.__class__.__name__))
1087 return '{}({})'.format(self.__class__.__name__, ', '.join(tokens))
1089 def __add__(self, other):
1090 atoms = self.copy()
1091 atoms += other
1092 return atoms
1094 def extend(self, other):
1095 """Extend atoms object by appending atoms from *other*."""
1096 if isinstance(other, Atom):
1097 other = self.__class__([other])
1099 n1 = len(self)
1100 n2 = len(other)
1102 for name, a1 in self.arrays.items():
1103 a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype)
1104 a[:n1] = a1
1105 if name == 'masses':
1106 a2 = other.get_masses()
1107 else:
1108 a2 = other.arrays.get(name)
1109 if a2 is not None:
1110 a[n1:] = a2
1111 self.arrays[name] = a
1113 for name, a2 in other.arrays.items():
1114 if name in self.arrays:
1115 continue
1116 a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype)
1117 a[n1:] = a2
1118 if name == 'masses':
1119 a[:n1] = self.get_masses()[:n1]
1120 else:
1121 a[:n1] = 0
1123 self.set_array(name, a)
1125 def __iadd__(self, other):
1126 self.extend(other)
1127 return self
1129 def append(self, atom):
1130 """Append atom to end."""
1131 self.extend(self.__class__([atom]))
1133 def __iter__(self):
1134 for i in range(len(self)):
1135 yield self[i]
1137 def __getitem__(self, i):
1138 """Return a subset of the atoms.
1140 i -- scalar integer, list of integers, or slice object
1141 describing which atoms to return.
1143 If i is a scalar, return an Atom object. If i is a list or a
1144 slice, return an Atoms object with the same cell, pbc, and
1145 other associated info as the original Atoms object. The
1146 indices of the constraints will be shuffled so that they match
1147 the indexing in the subset returned.
1149 """
1151 if isinstance(i, numbers.Integral):
1152 natoms = len(self)
1153 if i < -natoms or i >= natoms:
1154 raise IndexError('Index out of range.')
1156 return Atom(atoms=self, index=i)
1157 elif not isinstance(i, slice):
1158 i = np.array(i)
1159 if len(i) == 0:
1160 i = np.array([], dtype=int)
1161 # if i is a mask
1162 if i.dtype == bool:
1163 if len(i) != len(self):
1164 raise IndexError('Length of mask {} must equal '
1165 'number of atoms {}'
1166 .format(len(i), len(self)))
1167 i = np.arange(len(self))[i]
1169 import copy
1171 conadd = []
1172 # Constraints need to be deepcopied, but only the relevant ones.
1173 for con in copy.deepcopy(self.constraints):
1174 try:
1175 con.index_shuffle(self, i)
1176 except (IndexError, NotImplementedError):
1177 pass
1178 else:
1179 conadd.append(con)
1181 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1182 # should be communicated to the slice as well
1183 celldisp=self._celldisp)
1184 # TODO: Do we need to shuffle indices in adsorbate_info too?
1186 atoms.arrays = {}
1187 for name, a in self.arrays.items():
1188 atoms.arrays[name] = a[i].copy()
1190 atoms.constraints = conadd
1191 return atoms
1193 def __delitem__(self, i):
1194 from ase.constraints import FixAtoms
1195 for c in self._constraints:
1196 if not isinstance(c, FixAtoms):
1197 raise RuntimeError('Remove constraint using set_constraint() '
1198 'before deleting atoms.')
1200 if isinstance(i, list) and len(i) > 0:
1201 # Make sure a list of booleans will work correctly and not be
1202 # interpreted at 0 and 1 indices.
1203 i = np.array(i)
1205 if len(self._constraints) > 0:
1206 n = len(self)
1207 i = np.arange(n)[i]
1208 if isinstance(i, int):
1209 i = [i]
1210 constraints = []
1211 for c in self._constraints:
1212 c = c.delete_atoms(i, n)
1213 if c is not None:
1214 constraints.append(c)
1215 self.constraints = constraints
1217 mask = np.ones(len(self), bool)
1218 mask[i] = False
1219 for name, a in self.arrays.items():
1220 self.arrays[name] = a[mask]
1222 def pop(self, i=-1):
1223 """Remove and return atom at index *i* (default last)."""
1224 atom = self[i]
1225 atom.cut_reference_to_atoms()
1226 del self[i]
1227 return atom
1229 def __imul__(self, m):
1230 """In-place repeat of atoms."""
1231 if isinstance(m, int):
1232 m = (m, m, m)
1234 for x, vec in zip(m, self.cell):
1235 if x != 1 and not vec.any():
1236 raise ValueError('Cannot repeat along undefined lattice '
1237 'vector')
1239 M = np.prod(m)
1240 n = len(self)
1242 for name, a in self.arrays.items():
1243 self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1))
1245 positions = self.arrays['positions']
1246 i0 = 0
1247 for m0 in range(m[0]):
1248 for m1 in range(m[1]):
1249 for m2 in range(m[2]):
1250 i1 = i0 + n
1251 positions[i0:i1] += np.dot((m0, m1, m2), self.cell)
1252 i0 = i1
1254 if self.constraints is not None:
1255 self.constraints = [c.repeat(m, n) for c in self.constraints]
1257 self.cell = np.array([m[c] * self.cell[c] for c in range(3)])
1259 return self
1261 def repeat(self, rep):
1262 """Create new repeated atoms object.
1264 The *rep* argument should be a sequence of three positive
1265 integers like *(2,3,1)* or a single integer (*r*) equivalent
1266 to *(r,r,r)*."""
1268 atoms = self.copy()
1269 atoms *= rep
1270 return atoms
1272 def __mul__(self, rep):
1273 return self.repeat(rep)
1275 def translate(self, displacement):
1276 """Translate atomic positions.
1278 The displacement argument can be a float an xyz vector or an
1279 nx3 array (where n is the number of atoms)."""
1281 self.arrays['positions'] += np.array(displacement)
1283 def center(self, vacuum=None, axis=(0, 1, 2), about=None):
1284 """Center atoms in unit cell.
1286 Centers the atoms in the unit cell, so there is the same
1287 amount of vacuum on all sides.
1289 vacuum: float (default: None)
1290 If specified adjust the amount of vacuum when centering.
1291 If vacuum=10.0 there will thus be 10 Angstrom of vacuum
1292 on each side.
1293 axis: int or sequence of ints
1294 Axis or axes to act on. Default: Act on all axes.
1295 about: float or array (default: None)
1296 If specified, center the atoms about <about>.
1297 I.e., about=(0., 0., 0.) (or just "about=0.", interpreted
1298 identically), to center about the origin.
1299 """
1301 # Find the orientations of the faces of the unit cell
1302 cell = self.cell.complete()
1303 dirs = np.zeros_like(cell)
1305 lengths = cell.lengths()
1306 for i in range(3):
1307 dirs[i] = np.cross(cell[i - 1], cell[i - 2])
1308 dirs[i] /= np.linalg.norm(dirs[i])
1309 if dirs[i] @ cell[i] < 0.0:
1310 dirs[i] *= -1
1312 if isinstance(axis, int):
1313 axes = (axis,)
1314 else:
1315 axes = axis
1317 # Now, decide how much each basis vector should be made longer
1318 pos = self.positions
1319 longer = np.zeros(3)
1320 shift = np.zeros(3)
1321 for i in axes:
1322 if len(pos):
1323 scalarprod = pos @ dirs[i]
1324 p0 = scalarprod.min()
1325 p1 = scalarprod.max()
1326 else:
1327 p0 = 0
1328 p1 = 0
1329 height = cell[i] @ dirs[i]
1330 if vacuum is not None:
1331 lng = (p1 - p0 + 2 * vacuum) - height
1332 else:
1333 lng = 0.0 # Do not change unit cell size!
1334 top = lng + height - p1
1335 shf = 0.5 * (top - p0)
1336 cosphi = cell[i] @ dirs[i] / lengths[i]
1337 longer[i] = lng / cosphi
1338 shift[i] = shf / cosphi
1340 # Now, do it!
1341 translation = np.zeros(3)
1342 for i in axes:
1343 nowlen = lengths[i]
1344 if vacuum is not None:
1345 self.cell[i] = cell[i] * (1 + longer[i] / nowlen)
1346 translation += shift[i] * cell[i] / nowlen
1348 # We calculated translations using the completed cell,
1349 # so directions without cell vectors will have been centered
1350 # along a "fake" vector of length 1.
1351 # Therefore, we adjust by -0.5:
1352 if not any(self.cell[i]):
1353 translation[i] -= 0.5
1355 # Optionally, translate to center about a point in space.
1356 if about is not None:
1357 for n, vector in enumerate(self.cell):
1358 if n in axes:
1359 translation -= vector / 2.0
1360 translation[n] += about[n]
1362 self.positions += translation
1364 def get_center_of_mass(self, scaled=False, indices=None):
1365 """Get the center of mass.
1367 Parameters
1368 ----------
1369 scaled : bool
1370 If True, the center of mass in scaled coordinates is returned.
1371 indices : list | slice | str, default: None
1372 If specified, the center of mass of a subset of atoms is returned.
1373 """
1374 if indices is None:
1375 indices = slice(None)
1376 elif isinstance(indices, str):
1377 indices = string2index(indices)
1379 masses = self.get_masses()[indices]
1380 com = masses @ self.positions[indices] / masses.sum()
1381 if scaled:
1382 return self.cell.scaled_positions(com)
1383 return com # Cartesian coordinates
1385 def set_center_of_mass(self, com, scaled=False):
1386 """Set the center of mass.
1388 If scaled=True the center of mass is expected in scaled coordinates.
1389 Constraints are considered for scaled=False.
1390 """
1391 old_com = self.get_center_of_mass(scaled=scaled)
1392 difference = com - old_com
1393 if scaled:
1394 self.set_scaled_positions(self.get_scaled_positions() + difference)
1395 else:
1396 self.set_positions(self.get_positions() + difference)
1398 def get_moments_of_inertia(self, vectors=False):
1399 """Get the moments of inertia along the principal axes.
1401 The three principal moments of inertia are computed from the
1402 eigenvalues of the symmetric inertial tensor. Periodic boundary
1403 conditions are ignored. Units of the moments of inertia are
1404 amu*angstrom**2.
1405 """
1406 com = self.get_center_of_mass()
1407 positions = self.get_positions()
1408 positions -= com # translate center of mass to origin
1409 masses = self.get_masses()
1411 # Initialize elements of the inertial tensor
1412 I11 = I22 = I33 = I12 = I13 = I23 = 0.0
1413 for i in range(len(self)):
1414 x, y, z = positions[i]
1415 m = masses[i]
1417 I11 += m * (y ** 2 + z ** 2)
1418 I22 += m * (x ** 2 + z ** 2)
1419 I33 += m * (x ** 2 + y ** 2)
1420 I12 += -m * x * y
1421 I13 += -m * x * z
1422 I23 += -m * y * z
1424 Itensor = np.array([[I11, I12, I13],
1425 [I12, I22, I23],
1426 [I13, I23, I33]])
1428 evals, evecs = np.linalg.eigh(Itensor)
1429 if vectors:
1430 return evals, evecs.transpose()
1431 else:
1432 return evals
1434 def get_angular_momentum(self):
1435 """Get total angular momentum with respect to the center of mass."""
1436 com = self.get_center_of_mass()
1437 positions = self.get_positions()
1438 positions -= com # translate center of mass to origin
1439 return np.cross(positions, self.get_momenta()).sum(0)
1441 def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False):
1442 """Rotate atoms based on a vector and an angle, or two vectors.
1444 Parameters:
1446 a = None:
1447 Angle that the atoms is rotated around the vector 'v'. 'a'
1448 can also be a vector and then 'a' is rotated
1449 into 'v'.
1451 v:
1452 Vector to rotate the atoms around. Vectors can be given as
1453 strings: 'x', '-x', 'y', ... .
1455 center = (0, 0, 0):
1456 The center is kept fixed under the rotation. Use 'COM' to fix
1457 the center of mass, 'COP' to fix the center of positions or
1458 'COU' to fix the center of cell.
1460 rotate_cell = False:
1461 If true the cell is also rotated.
1463 Examples:
1465 Rotate 90 degrees around the z-axis, so that the x-axis is
1466 rotated into the y-axis:
1468 >>> atoms = Atoms()
1469 >>> atoms.rotate(90, 'z')
1470 >>> atoms.rotate(90, (0, 0, 1))
1471 >>> atoms.rotate(-90, '-z')
1472 >>> atoms.rotate('x', 'y')
1473 >>> atoms.rotate((1, 0, 0), (0, 1, 0))
1474 """
1476 if not isinstance(a, numbers.Real):
1477 a, v = v, a
1479 norm = np.linalg.norm
1480 v = string2vector(v)
1482 normv = norm(v)
1484 if normv == 0.0:
1485 raise ZeroDivisionError('Cannot rotate: norm(v) == 0')
1487 if isinstance(a, numbers.Real):
1488 a *= pi / 180
1489 v /= normv
1490 c = cos(a)
1491 s = sin(a)
1492 else:
1493 v2 = string2vector(a)
1494 v /= normv
1495 normv2 = np.linalg.norm(v2)
1496 if normv2 == 0:
1497 raise ZeroDivisionError('Cannot rotate: norm(a) == 0')
1498 v2 /= norm(v2)
1499 c = np.dot(v, v2)
1500 v = np.cross(v, v2)
1501 s = norm(v)
1502 # In case *v* and *a* are parallel, np.cross(v, v2) vanish
1503 # and can't be used as a rotation axis. However, in this
1504 # case any rotation axis perpendicular to v2 will do.
1505 eps = 1e-7
1506 if s < eps:
1507 v = np.cross((0, 0, 1), v2)
1508 if norm(v) < eps:
1509 v = np.cross((1, 0, 0), v2)
1510 assert norm(v) >= eps
1511 elif s > 0:
1512 v /= s
1514 center = self._centering_as_array(center)
1516 p = self.arrays['positions'] - center
1517 self.arrays['positions'][:] = (c * p -
1518 np.cross(p, s * v) +
1519 np.outer(np.dot(p, v), (1.0 - c) * v) +
1520 center)
1521 if rotate_cell:
1522 rotcell = self.get_cell()
1523 rotcell[:] = (c * rotcell -
1524 np.cross(rotcell, s * v) +
1525 np.outer(np.dot(rotcell, v), (1.0 - c) * v))
1526 self.set_cell(rotcell)
1528 def _centering_as_array(self, center):
1529 if isinstance(center, str):
1530 if center.lower() == 'com':
1531 center = self.get_center_of_mass()
1532 elif center.lower() == 'cop':
1533 center = self.get_positions().mean(axis=0)
1534 elif center.lower() == 'cou':
1535 center = self.get_cell().sum(axis=0) / 2
1536 else:
1537 raise ValueError('Cannot interpret center')
1538 else:
1539 center = np.array(center, float)
1540 return center
1542 def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)):
1543 """Rotate atoms via Euler angles (in degrees).
1545 See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation.
1547 Parameters:
1549 center :
1550 The point to rotate about. A sequence of length 3 with the
1551 coordinates, or 'COM' to select the center of mass, 'COP' to
1552 select center of positions or 'COU' to select center of cell.
1553 phi :
1554 The 1st rotation angle around the z axis.
1555 theta :
1556 Rotation around the x axis.
1557 psi :
1558 2nd rotation around the z axis.
1560 """
1561 center = self._centering_as_array(center)
1563 phi *= pi / 180
1564 theta *= pi / 180
1565 psi *= pi / 180
1567 # First move the molecule to the origin In contrast to MATLAB,
1568 # numpy broadcasts the smaller array to the larger row-wise,
1569 # so there is no need to play with the Kronecker product.
1570 rcoords = self.positions - center
1571 # First Euler rotation about z in matrix form
1572 D = np.array(((cos(phi), sin(phi), 0.),
1573 (-sin(phi), cos(phi), 0.),
1574 (0., 0., 1.)))
1575 # Second Euler rotation about x:
1576 C = np.array(((1., 0., 0.),
1577 (0., cos(theta), sin(theta)),
1578 (0., -sin(theta), cos(theta))))
1579 # Third Euler rotation, 2nd rotation about z:
1580 B = np.array(((cos(psi), sin(psi), 0.),
1581 (-sin(psi), cos(psi), 0.),
1582 (0., 0., 1.)))
1583 # Total Euler rotation
1584 A = np.dot(B, np.dot(C, D))
1585 # Do the rotation
1586 rcoords = np.dot(A, np.transpose(rcoords))
1587 # Move back to the rotation point
1588 self.positions = np.transpose(rcoords) + center
1590 def get_dihedral(self, a0, a1, a2, a3, mic=False):
1591 """Calculate dihedral angle.
1593 Calculate dihedral angle (in degrees) between the vectors a0->a1
1594 and a2->a3.
1596 Use mic=True to use the Minimum Image Convention and calculate the
1597 angle across periodic boundaries.
1598 """
1599 return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0]
1601 def get_dihedrals(self, indices, mic=False):
1602 """Calculate dihedral angles.
1604 Calculate dihedral angles (in degrees) between the list of vectors
1605 a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices.
1607 Use mic=True to use the Minimum Image Convention and calculate the
1608 angles across periodic boundaries.
1609 """
1610 from ase.geometry import get_dihedrals
1612 indices = np.array(indices)
1613 assert indices.shape[1] == 4
1615 a0s = self.positions[indices[:, 0]]
1616 a1s = self.positions[indices[:, 1]]
1617 a2s = self.positions[indices[:, 2]]
1618 a3s = self.positions[indices[:, 3]]
1620 # vectors 0->1, 1->2, 2->3
1621 v0 = a1s - a0s
1622 v1 = a2s - a1s
1623 v2 = a3s - a2s
1625 cell = None
1626 pbc = None
1628 if mic:
1629 cell = self.cell
1630 pbc = self.pbc
1632 return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc)
1634 def _masked_rotate(self, center, axis, diff, mask):
1635 # do rotation of subgroup by copying it to temporary atoms object
1636 # and then rotating that
1637 #
1638 # recursive object definition might not be the most elegant thing,
1639 # more generally useful might be a rotation function with a mask?
1640 group = self.__class__()
1641 for i in range(len(self)):
1642 if mask[i]:
1643 group += self[i]
1644 group.translate(-center)
1645 group.rotate(diff * 180 / pi, axis)
1646 group.translate(center)
1647 # set positions in original atoms object
1648 j = 0
1649 for i in range(len(self)):
1650 if mask[i]:
1651 self.positions[i] = group[j].position
1652 j += 1
1654 def set_dihedral(self, a1, a2, a3, a4, angle,
1655 mask=None, indices=None):
1656 """Set the dihedral angle (degrees) between vectors a1->a2 and
1657 a3->a4 by changing the atom indexed by a4.
1659 If mask is not None, all the atoms described in mask
1660 (read: the entire subgroup) are moved. Alternatively to the mask,
1661 the indices of the atoms to be rotated can be supplied. If both
1662 *mask* and *indices* are given, *indices* overwrites *mask*.
1664 **Important**: If *mask* or *indices* is given and does not contain
1665 *a4*, *a4* will NOT be moved. In most cases you therefore want
1666 to include *a4* in *mask*/*indices*.
1668 Example: the following defines a very crude
1669 ethane-like molecule and twists one half of it by 30 degrees.
1671 >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0],
1672 ... [1, 0, 0], [2, 1, 0], [2, -1, 0]])
1673 >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1])
1674 """
1676 angle *= pi / 180
1678 # if not provided, set mask to the last atom in the
1679 # dihedral description
1680 if mask is None and indices is None:
1681 mask = np.zeros(len(self))
1682 mask[a4] = 1
1683 elif indices is not None:
1684 mask = [index in indices for index in range(len(self))]
1686 # compute necessary in dihedral change, from current value
1687 current = self.get_dihedral(a1, a2, a3, a4) * pi / 180
1688 diff = angle - current
1689 axis = self.positions[a3] - self.positions[a2]
1690 center = self.positions[a3]
1691 self._masked_rotate(center, axis, diff, mask)
1693 def rotate_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None):
1694 """Rotate dihedral angle.
1696 Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a
1697 predefined dihedral angle, starting from its current configuration.
1698 """
1699 start = self.get_dihedral(a1, a2, a3, a4)
1700 self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices)
1702 def get_angle(self, a1, a2, a3, mic=False):
1703 """Get angle formed by three atoms.
1705 Calculate angle in degrees between the vectors a2->a1 and
1706 a2->a3.
1708 Use mic=True to use the Minimum Image Convention and calculate the
1709 angle across periodic boundaries.
1710 """
1711 return self.get_angles([[a1, a2, a3]], mic=mic)[0]
1713 def get_angles(self, indices, mic=False):
1714 """Get angle formed by three atoms for multiple groupings.
1716 Calculate angle in degrees between vectors between atoms a2->a1
1717 and a2->a3, where a1, a2, and a3 are in each row of indices.
1719 Use mic=True to use the Minimum Image Convention and calculate
1720 the angle across periodic boundaries.
1721 """
1722 from ase.geometry import get_angles
1724 indices = np.array(indices)
1725 assert indices.shape[1] == 3
1727 a1s = self.positions[indices[:, 0]]
1728 a2s = self.positions[indices[:, 1]]
1729 a3s = self.positions[indices[:, 2]]
1731 v12 = a1s - a2s
1732 v32 = a3s - a2s
1734 cell = None
1735 pbc = None
1737 if mic:
1738 cell = self.cell
1739 pbc = self.pbc
1741 return get_angles(v12, v32, cell=cell, pbc=pbc)
1743 def set_angle(self, a1, a2=None, a3=None, angle=None, mask=None,
1744 indices=None, add=False):
1745 """Set angle (in degrees) formed by three atoms.
1747 Sets the angle between vectors *a2*->*a1* and *a2*->*a3*.
1749 If *add* is `True`, the angle will be changed by the value given.
1751 Same usage as in :meth:`ase.Atoms.set_dihedral`.
1752 If *mask* and *indices*
1753 are given, *indices* overwrites *mask*. If *mask* and *indices*
1754 are not set, only *a3* is moved."""
1756 if any(a is None for a in [a2, a3, angle]):
1757 raise ValueError('a2, a3, and angle must not be None')
1759 # If not provided, set mask to the last atom in the angle description
1760 if mask is None and indices is None:
1761 mask = np.zeros(len(self))
1762 mask[a3] = 1
1763 elif indices is not None:
1764 mask = [index in indices for index in range(len(self))]
1766 if add:
1767 diff = angle
1768 else:
1769 # Compute necessary in angle change, from current value
1770 diff = angle - self.get_angle(a1, a2, a3)
1772 diff *= pi / 180
1773 # Do rotation of subgroup by copying it to temporary atoms object and
1774 # then rotating that
1775 v10 = self.positions[a1] - self.positions[a2]
1776 v12 = self.positions[a3] - self.positions[a2]
1777 v10 /= np.linalg.norm(v10)
1778 v12 /= np.linalg.norm(v12)
1779 axis = np.cross(v10, v12)
1780 center = self.positions[a2]
1781 self._masked_rotate(center, axis, diff, mask)
1783 def rattle(self, stdev=0.001, seed=None, rng=None):
1784 """Randomly displace atoms.
1786 This method adds random displacements to the atomic positions,
1787 taking a possible constraint into account. The random numbers are
1788 drawn from a normal distribution of standard deviation stdev.
1790 By default, the random number generator always uses the same seed (42)
1791 for repeatability. You can provide your own seed (an integer), or if you
1792 want the randomness to be different each time you run a script, then
1793 provide `rng=numpy.random`. For a parallel calculation, it is important
1794 to use the same seed on all processors! """
1796 if seed is not None and rng is not None:
1797 raise ValueError('Please do not provide both seed and rng.')
1799 if rng is None:
1800 if seed is None:
1801 seed = 42
1802 rng = np.random.RandomState(seed)
1803 positions = self.arrays['positions']
1804 self.set_positions(positions +
1805 rng.normal(scale=stdev, size=positions.shape))
1807 def get_distance(self, a0, a1, mic=False, vector=False):
1808 """Return distance between two atoms.
1810 Use mic=True to use the Minimum Image Convention.
1811 vector=True gives the distance vector (from a0 to a1).
1812 """
1813 return self.get_distances(a0, [a1], mic=mic, vector=vector)[0]
1815 def get_distances(self, a, indices, mic=False, vector=False):
1816 """Return distances of atom No.i with a list of atoms.
1818 Use mic=True to use the Minimum Image Convention.
1819 vector=True gives the distance vector (from a to self[indices]).
1820 """
1821 from ase.geometry import get_distances
1823 R = self.arrays['positions']
1824 p1 = [R[a]]
1825 p2 = R[indices]
1827 cell = None
1828 pbc = None
1830 if mic:
1831 cell = self.cell
1832 pbc = self.pbc
1834 D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc)
1836 if vector:
1837 D.shape = (-1, 3)
1838 return D
1839 else:
1840 D_len.shape = (-1,)
1841 return D_len
1843 def get_all_distances(self, mic=False, vector=False):
1844 """Return distances of all of the atoms with all of the atoms.
1846 Use mic=True to use the Minimum Image Convention.
1847 """
1848 from ase.geometry import get_distances
1850 R = self.arrays['positions']
1852 cell = None
1853 pbc = None
1855 if mic:
1856 cell = self.cell
1857 pbc = self.pbc
1859 D, D_len = get_distances(R, cell=cell, pbc=pbc)
1861 if vector:
1862 return D
1863 else:
1864 return D_len
1866 def set_distance(self, a0, a1, distance, fix=0.5, mic=False,
1867 mask=None, indices=None, add=False, factor=False):
1868 """Set the distance between two atoms.
1870 Set the distance between atoms *a0* and *a1* to *distance*.
1871 By default, the center of the two atoms will be fixed. Use
1872 *fix=0* to fix the first atom, *fix=1* to fix the second
1873 atom and *fix=0.5* (default) to fix the center of the bond.
1875 If *mask* or *indices* are set (*mask* overwrites *indices*),
1876 only the atoms defined there are moved
1877 (see :meth:`ase.Atoms.set_dihedral`).
1879 When *add* is true, the distance is changed by the value given.
1880 In combination
1881 with *factor* True, the value given is a factor scaling the distance.
1883 It is assumed that the atoms in *mask*/*indices* move together
1884 with *a1*. If *fix=1*, only *a0* will therefore be moved."""
1885 from ase.geometry import find_mic
1887 if a0 % len(self) == a1 % len(self):
1888 raise ValueError('a0 and a1 must not be the same')
1890 if add:
1891 oldDist = self.get_distance(a0, a1, mic=mic)
1892 if factor:
1893 newDist = oldDist * distance
1894 else:
1895 newDist = oldDist + distance
1896 self.set_distance(a0, a1, newDist, fix=fix, mic=mic,
1897 mask=mask, indices=indices, add=False,
1898 factor=False)
1899 return
1901 R = self.arrays['positions']
1902 D = np.array([R[a1] - R[a0]])
1904 if mic:
1905 D, D_len = find_mic(D, self.cell, self.pbc)
1906 else:
1907 D_len = np.array([np.sqrt((D**2).sum())])
1908 x = 1.0 - distance / D_len[0]
1910 if mask is None and indices is None:
1911 indices = [a0, a1]
1912 elif mask:
1913 indices = [i for i in range(len(self)) if mask[i]]
1915 for i in indices:
1916 if i == a0:
1917 R[a0] += (x * fix) * D[0]
1918 else:
1919 R[i] -= (x * (1.0 - fix)) * D[0]
1921 def get_scaled_positions(self, wrap=True):
1922 """Get positions relative to unit cell.
1924 If wrap is True, atoms outside the unit cell will be wrapped into
1925 the cell in those directions with periodic boundary conditions
1926 so that the scaled coordinates are between zero and one.
1928 If any cell vectors are zero, the corresponding coordinates
1929 are evaluated as if the cell were completed using
1930 ``cell.complete()``. This means coordinates will be Cartesian
1931 as long as the non-zero cell vectors span a Cartesian axis or
1932 plane."""
1934 fractional = self.cell.scaled_positions(self.positions)
1936 if wrap:
1937 for i, periodic in enumerate(self.pbc):
1938 if periodic:
1939 # Yes, we need to do it twice.
1940 # See the scaled_positions.py test.
1941 fractional[:, i] %= 1.0
1942 fractional[:, i] %= 1.0
1944 return fractional
1946 def set_scaled_positions(self, scaled):
1947 """Set positions relative to unit cell."""
1948 self.positions[:] = self.cell.cartesian_positions(scaled)
1950 def wrap(self, **wrap_kw):
1951 """Wrap positions to unit cell.
1953 Parameters:
1955 wrap_kw: (keyword=value) pairs
1956 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
1957 see :func:`ase.geometry.wrap_positions`
1958 """
1960 if 'pbc' not in wrap_kw:
1961 wrap_kw['pbc'] = self.pbc
1963 self.positions[:] = self.get_positions(wrap=True, **wrap_kw)
1965 def get_temperature(self):
1966 """Get the temperature in Kelvin."""
1967 ekin = self.get_kinetic_energy()
1968 return 2 * ekin / (self.get_number_of_degrees_of_freedom() * units.kB)
1970 def __eq__(self, other):
1971 """Check for identity of two atoms objects.
1973 Identity means: same positions, atomic numbers, unit cell and
1974 periodic boundary conditions."""
1975 if not isinstance(other, Atoms):
1976 return False
1977 a = self.arrays
1978 b = other.arrays
1979 return (len(self) == len(other) and
1980 (a['positions'] == b['positions']).all() and
1981 (a['numbers'] == b['numbers']).all() and
1982 (self.cell == other.cell).all() and
1983 (self.pbc == other.pbc).all())
1985 def __ne__(self, other):
1986 """Check if two atoms objects are not equal.
1988 Any differences in positions, atomic numbers, unit cell or
1989 periodic boundary condtions make atoms objects not equal.
1990 """
1991 eq = self.__eq__(other)
1992 if eq is NotImplemented:
1993 return eq
1994 else:
1995 return not eq
1997 # @deprecated('Please use atoms.cell.volume')
1998 # We kind of want to deprecate this, but the ValueError behaviour
1999 # might be desirable. Should we do this?
2000 def get_volume(self):
2001 """Get volume of unit cell."""
2002 if self.cell.rank != 3:
2003 raise ValueError(
2004 'You have {} lattice vectors: volume not defined'
2005 .format(self.cell.rank))
2006 return self.cell.volume
2008 def _get_positions(self):
2009 """Return reference to positions-array for in-place manipulations."""
2010 return self.arrays['positions']
2012 def _set_positions(self, pos):
2013 """Set positions directly, bypassing constraints."""
2014 self.arrays['positions'][:] = pos
2016 positions = property(_get_positions, _set_positions,
2017 doc='Attribute for direct ' +
2018 'manipulation of the positions.')
2020 def _get_atomic_numbers(self):
2021 """Return reference to atomic numbers for in-place
2022 manipulations."""
2023 return self.arrays['numbers']
2025 numbers = property(_get_atomic_numbers, set_atomic_numbers,
2026 doc='Attribute for direct ' +
2027 'manipulation of the atomic numbers.')
2029 @property
2030 def cell(self):
2031 """The :class:`ase.cell.Cell` for direct manipulation."""
2032 return self._cellobj
2034 @cell.setter
2035 def cell(self, cell):
2036 cell = Cell.ascell(cell)
2037 self._cellobj[:] = cell
2039 def write(self, filename, format=None, **kwargs):
2040 """Write atoms object to a file.
2042 see ase.io.write for formats.
2043 kwargs are passed to ase.io.write.
2044 """
2045 from ase.io import write
2046 write(filename, self, format, **kwargs)
2048 def iterimages(self):
2049 yield self
2051 def __ase_optimizable__(self):
2052 from ase.optimize.optimize import OptimizableAtoms
2053 return OptimizableAtoms(self)
2055 def edit(self):
2056 """Modify atoms interactively through ASE's GUI viewer.
2058 Conflicts leading to undesirable behaviour might arise
2059 when matplotlib has been pre-imported with certain
2060 incompatible backends and while trying to use the
2061 plot feature inside the interactive GUI. To circumvent,
2062 please set matplotlib.use('gtk') before calling this
2063 method.
2064 """
2065 from ase.gui.gui import GUI
2066 from ase.gui.images import Images
2067 images = Images([self])
2068 gui = GUI(images)
2069 gui.run()
2072def string2vector(v):
2073 if isinstance(v, str):
2074 if v[0] == '-':
2075 return -string2vector(v[1:])
2076 w = np.zeros(3)
2077 w['xyz'.index(v)] = 1.0
2078 return w
2079 return np.array(v, float)
2082def default(data, dflt):
2083 """Helper function for setting default values."""
2084 if data is None:
2085 return None
2086 elif isinstance(data, (list, tuple)):
2087 newdata = []
2088 allnone = True
2089 for x in data:
2090 if x is None:
2091 newdata.append(dflt)
2092 else:
2093 newdata.append(x)
2094 allnone = False
2095 if allnone:
2096 return None
2097 return newdata
2098 else:
2099 return data