Coverage for /builds/hweiske/ase/ase/atoms.py: 93.89%
966 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +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:
40 symbols: str (formula) or list of str
41 Can be a string formula, a list of symbols or a list of
42 Atom objects. Examples: 'H2O', 'COPt12', ['H', 'H', 'O'],
43 [Atom('Ne', (x, y, z)), ...].
44 positions: list of xyz-positions
45 Atomic positions. Anything that can be converted to an
46 ndarray of shape (n, 3) will do: [(x1,y1,z1), (x2,y2,z2),
47 ...].
48 scaled_positions: list of scaled-positions
49 Like positions, but given in units of the unit cell.
50 Can not be set at the same time as positions.
51 numbers: list of int
52 Atomic numbers (use only one of symbols/numbers).
53 tags: list of int
54 Special purpose tags.
55 momenta: list of xyz-momenta
56 Momenta for all atoms.
57 masses: list of float
58 Atomic masses in atomic units.
59 magmoms: list of float or list of xyz-values
60 Magnetic moments. Can be either a single value for each atom
61 for collinear calculations or three numbers for each atom for
62 non-collinear calculations.
63 charges: list of float
64 Initial atomic charges.
65 cell: 3x3 matrix or length 3 or 6 vector
66 Unit cell vectors. Can also be given as just three
67 numbers for orthorhombic cells, or 6 numbers, where
68 first three are lengths of unit cell vectors, and the
69 other three are angles between them (in degrees), in following order:
70 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)].
71 First vector will lie in x-direction, second in xy-plane,
72 and the third one in z-positive subspace.
73 Default value: [0, 0, 0].
74 celldisp: Vector
75 Unit cell displacement vector. To visualize a displaced cell
76 around the center of mass of a Systems of atoms. Default value
77 = (0,0,0)
78 pbc: one or three bool
79 Periodic boundary conditions flags. Examples: True,
80 False, 0, 1, (1, 1, 0), (True, False, False). Default
81 value: False.
82 constraint: constraint object(s)
83 Used for applying one or more constraints during structure
84 optimization.
85 calculator: calculator object
86 Used to attach a calculator for calculating energies and atomic
87 forces.
88 info: dict of key-value pairs
89 Dictionary of key-value pairs with additional information
90 about the system. The following keys may be used by ase:
92 - spacegroup: Spacegroup instance
93 - unit_cell: 'conventional' | 'primitive' | int | 3 ints
94 - adsorbate_info: Information about special adsorption sites
96 Items in the info attribute survives copy and slicing and can
97 be stored in and retrieved from trajectory files given that the
98 key is a string, the value is JSON-compatible and, if the value is a
99 user-defined object, its base class is importable. One should
100 not make any assumptions about the existence of keys.
102 Examples:
104 These three are equivalent:
106 >>> from ase import Atom
108 >>> d = 1.104 # N2 bondlength
109 >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)])
110 >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)])
111 >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))])
113 FCC gold:
115 >>> a = 4.05 # Gold lattice constant
116 >>> b = a / 2
117 >>> fcc = Atoms('Au',
118 ... cell=[(0, b, b), (b, 0, b), (b, b, 0)],
119 ... pbc=True)
121 Hydrogen wire:
123 >>> d = 0.9 # H-H distance
124 >>> h = Atoms('H', positions=[(0, 0, 0)],
125 ... cell=(d, 0, 0),
126 ... pbc=(1, 0, 0))
127 """
129 ase_objtype = 'atoms' # For JSONability
131 def __init__(self, symbols=None,
132 positions=None, numbers=None,
133 tags=None, momenta=None, masses=None,
134 magmoms=None, charges=None,
135 scaled_positions=None,
136 cell=None, pbc=None, celldisp=None,
137 constraint=None,
138 calculator=None,
139 info=None,
140 velocities=None):
142 self._cellobj = Cell.new()
143 self._pbc = np.zeros(3, bool)
145 atoms = None
147 if hasattr(symbols, 'get_positions'):
148 atoms = symbols
149 symbols = None
150 elif (isinstance(symbols, (list, tuple)) and
151 len(symbols) > 0 and isinstance(symbols[0], Atom)):
152 # Get data from a list or tuple of Atom objects:
153 data = [[atom.get_raw(name) for atom in symbols]
154 for name in
155 ['position', 'number', 'tag', 'momentum',
156 'mass', 'magmom', 'charge']]
157 atoms = self.__class__(None, *data)
158 symbols = None
160 if atoms is not None:
161 # Get data from another Atoms object:
162 if scaled_positions is not None:
163 raise NotImplementedError
164 if symbols is None and numbers is None:
165 numbers = atoms.get_atomic_numbers()
166 if positions is None:
167 positions = atoms.get_positions()
168 if tags is None and atoms.has('tags'):
169 tags = atoms.get_tags()
170 if momenta is None and atoms.has('momenta'):
171 momenta = atoms.get_momenta()
172 if magmoms is None and atoms.has('initial_magmoms'):
173 magmoms = atoms.get_initial_magnetic_moments()
174 if masses is None and atoms.has('masses'):
175 masses = atoms.get_masses()
176 if charges is None and atoms.has('initial_charges'):
177 charges = atoms.get_initial_charges()
178 if cell is None:
179 cell = atoms.get_cell()
180 if celldisp is None:
181 celldisp = atoms.get_celldisp()
182 if pbc is None:
183 pbc = atoms.get_pbc()
184 if constraint is None:
185 constraint = [c.copy() for c in atoms.constraints]
186 if calculator is None:
187 calculator = atoms.calc
188 if info is None:
189 info = copy.deepcopy(atoms.info)
191 self.arrays = {}
193 if symbols is None:
194 if numbers is None:
195 if positions is not None:
196 natoms = len(positions)
197 elif scaled_positions is not None:
198 natoms = len(scaled_positions)
199 else:
200 natoms = 0
201 numbers = np.zeros(natoms, int)
202 self.new_array('numbers', numbers, int)
203 else:
204 if numbers is not None:
205 raise TypeError(
206 'Use only one of "symbols" and "numbers".')
207 else:
208 self.new_array('numbers', symbols2numbers(symbols), int)
210 if self.numbers.ndim != 1:
211 raise ValueError('"numbers" must be 1-dimensional.')
213 if cell is None:
214 cell = np.zeros((3, 3))
215 self.set_cell(cell)
217 if celldisp is None:
218 celldisp = np.zeros(shape=(3, 1))
219 self.set_celldisp(celldisp)
221 if positions is None:
222 if scaled_positions is None:
223 positions = np.zeros((len(self.arrays['numbers']), 3))
224 else:
225 assert self.cell.rank == 3
226 positions = np.dot(scaled_positions, self.cell)
227 else:
228 if scaled_positions is not None:
229 raise TypeError(
230 'Use only one of "symbols" and "numbers".')
231 self.new_array('positions', positions, float, (3,))
233 self.set_constraint(constraint)
234 self.set_tags(default(tags, 0))
235 self.set_masses(default(masses, None))
236 self.set_initial_magnetic_moments(default(magmoms, 0.0))
237 self.set_initial_charges(default(charges, 0.0))
238 if pbc is None:
239 pbc = False
240 self.set_pbc(pbc)
241 self.set_momenta(default(momenta, (0.0, 0.0, 0.0)),
242 apply_constraint=False)
244 if velocities is not None:
245 if momenta is None:
246 self.set_velocities(velocities)
247 else:
248 raise TypeError(
249 'Use only one of "momenta" and "velocities".')
251 if info is None:
252 self.info = {}
253 else:
254 self.info = dict(info)
256 self.calc = calculator
258 @property
259 def symbols(self):
260 """Get chemical symbols as a :class:`ase.symbols.Symbols` object.
262 The object works like ``atoms.numbers`` except its values
263 are strings. It supports in-place editing."""
264 return Symbols(self.numbers)
266 @symbols.setter
267 def symbols(self, obj):
268 new_symbols = Symbols.fromsymbols(obj)
269 self.numbers[:] = new_symbols.numbers
271 @deprecated("Please use atoms.calc = calc", DeprecationWarning)
272 def set_calculator(self, calc=None):
273 """Attach calculator object.
275 .. deprecated:: 3.20.0
276 Please use the equivalent ``atoms.calc = calc`` instead of this
277 method.
278 """
280 self.calc = calc
282 @deprecated("Please use atoms.calc", DeprecationWarning)
283 def get_calculator(self):
284 """Get currently attached calculator object.
286 .. deprecated:: 3.20.0
287 Please use the equivalent ``atoms.calc`` instead of
288 ``atoms.get_calculator()``.
289 """
291 return self.calc
293 @property
294 def calc(self):
295 """Calculator object."""
296 return self._calc
298 @calc.setter
299 def calc(self, calc):
300 self._calc = calc
301 if hasattr(calc, 'set_atoms'):
302 calc.set_atoms(self)
304 @calc.deleter
305 @deprecated('Please use atoms.calc = None', DeprecationWarning)
306 def calc(self):
307 """Delete calculator
309 .. deprecated:: 3.20.0
310 Please use ``atoms.calc = None``
311 """
312 self._calc = None
314 @property
315 @deprecated('Please use atoms.cell.rank instead', DeprecationWarning)
316 def number_of_lattice_vectors(self):
317 """Number of (non-zero) lattice vectors.
319 .. deprecated:: 3.21.0
320 Please use ``atoms.cell.rank`` instead
321 """
322 return self.cell.rank
324 def set_constraint(self, constraint=None):
325 """Apply one or more constrains.
327 The *constraint* argument must be one constraint object or a
328 list of constraint objects."""
329 if constraint is None:
330 self._constraints = []
331 else:
332 if isinstance(constraint, list):
333 self._constraints = constraint
334 elif isinstance(constraint, tuple):
335 self._constraints = list(constraint)
336 else:
337 self._constraints = [constraint]
339 def _get_constraints(self):
340 return self._constraints
342 def _del_constraints(self):
343 self._constraints = []
345 constraints = property(_get_constraints, set_constraint, _del_constraints,
346 'Constraints of the atoms.')
348 def set_cell(self, cell, scale_atoms=False, apply_constraint=True):
349 """Set unit cell vectors.
351 Parameters:
353 cell: 3x3 matrix or length 3 or 6 vector
354 Unit cell. A 3x3 matrix (the three unit cell vectors) or
355 just three numbers for an orthorhombic cell. Another option is
356 6 numbers, which describes unit cell with lengths of unit cell
357 vectors and with angles between them (in degrees), in following
358 order: [len(a), len(b), len(c), angle(b,c), angle(a,c),
359 angle(a,b)]. First vector will lie in x-direction, second in
360 xy-plane, and the third one in z-positive subspace.
361 scale_atoms: bool
362 Fix atomic positions or move atoms with the unit cell?
363 Default behavior is to *not* move the atoms (scale_atoms=False).
364 apply_constraint: bool
365 Whether to apply constraints to the given cell.
367 Examples:
369 Two equivalent ways to define an orthorhombic cell:
371 >>> atoms = Atoms('He')
372 >>> a, b, c = 7, 7.5, 8
373 >>> atoms.set_cell([a, b, c])
374 >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)])
376 FCC unit cell:
378 >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)])
380 Hexagonal unit cell:
382 >>> atoms.set_cell([a, a, c, 90, 90, 120])
384 Rhombohedral unit cell:
386 >>> alpha = 77
387 >>> atoms.set_cell([a, a, a, alpha, alpha, alpha])
388 """
390 # Override pbcs if and only if given a Cell object:
391 cell = Cell.new(cell)
393 # XXX not working well during initialize due to missing _constraints
394 if apply_constraint and hasattr(self, '_constraints'):
395 for constraint in self.constraints:
396 if hasattr(constraint, 'adjust_cell'):
397 constraint.adjust_cell(self, cell)
399 if scale_atoms:
400 M = np.linalg.solve(self.cell.complete(), cell.complete())
401 self.positions[:] = np.dot(self.positions, M)
403 self.cell[:] = cell
405 def set_celldisp(self, celldisp):
406 """Set the unit cell displacement vectors."""
407 celldisp = np.array(celldisp, float)
408 self._celldisp = celldisp
410 def get_celldisp(self):
411 """Get the unit cell displacement vectors."""
412 return self._celldisp.copy()
414 def get_cell(self, complete=False):
415 """Get the three unit cell vectors as a `class`:ase.cell.Cell` object.
417 The Cell object resembles a 3x3 ndarray, and cell[i, j]
418 is the jth Cartesian coordinate of the ith cell vector."""
419 if complete:
420 cell = self.cell.complete()
421 else:
422 cell = self.cell.copy()
424 return cell
426 @deprecated('Please use atoms.cell.cellpar() instead', DeprecationWarning)
427 def get_cell_lengths_and_angles(self):
428 """Get unit cell parameters. Sequence of 6 numbers.
430 First three are unit cell vector lengths and second three
431 are angles between them::
433 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
435 in degrees.
437 .. deprecated:: 3.21.0
438 Please use ``atoms.cell.cellpar()`` instead
439 """
440 return self.cell.cellpar()
442 @deprecated('Please use atoms.cell.reciprocal()', DeprecationWarning)
443 def get_reciprocal_cell(self):
444 """Get the three reciprocal lattice vectors as a 3x3 ndarray.
446 Note that the commonly used factor of 2 pi for Fourier
447 transforms is not included here.
449 .. deprecated:: 3.21.0
450 Please use ``atoms.cell.reciprocal()``
451 """
452 return self.cell.reciprocal()
454 @property
455 def pbc(self):
456 """Reference to pbc-flags for in-place manipulations."""
457 return self._pbc
459 @pbc.setter
460 def pbc(self, pbc):
461 self._pbc[:] = pbc
463 def set_pbc(self, pbc):
464 """Set periodic boundary condition flags."""
465 self.pbc = pbc
467 def get_pbc(self):
468 """Get periodic boundary condition flags."""
469 return self.pbc.copy()
471 def new_array(self, name, a, dtype=None, shape=None):
472 """Add new array.
474 If *shape* is not *None*, the shape of *a* will be checked."""
476 if dtype is not None:
477 a = np.array(a, dtype, order='C')
478 if len(a) == 0 and shape is not None:
479 a.shape = (-1,) + shape
480 else:
481 if not a.flags['C_CONTIGUOUS']:
482 a = np.ascontiguousarray(a)
483 else:
484 a = a.copy()
486 if name in self.arrays:
487 raise RuntimeError(f'Array {name} already present')
489 for b in self.arrays.values():
490 if len(a) != len(b):
491 raise ValueError('Array "%s" has wrong length: %d != %d.' %
492 (name, len(a), len(b)))
493 break
495 if shape is not None and a.shape[1:] != shape:
496 raise ValueError(
497 f'Array "{name}" has wrong shape {a.shape} != '
498 f'{(a.shape[0:1] + shape)}.')
500 self.arrays[name] = a
502 def get_array(self, name, copy=True):
503 """Get an array.
505 Returns a copy unless the optional argument copy is false.
506 """
507 if copy:
508 return self.arrays[name].copy()
509 else:
510 return self.arrays[name]
512 def set_array(self, name, a, dtype=None, shape=None):
513 """Update array.
515 If *shape* is not *None*, the shape of *a* will be checked.
516 If *a* is *None*, then the array is deleted."""
518 b = self.arrays.get(name)
519 if b is None:
520 if a is not None:
521 self.new_array(name, a, dtype, shape)
522 else:
523 if a is None:
524 del self.arrays[name]
525 else:
526 a = np.asarray(a)
527 if a.shape != b.shape:
528 raise ValueError(
529 f'Array "{name}" has wrong shape '
530 f'{a.shape} != {b.shape}.')
531 b[:] = a
533 def has(self, name):
534 """Check for existence of array.
536 name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms',
537 'initial_charges'."""
538 # XXX extend has to calculator properties
539 return name in self.arrays
541 def set_atomic_numbers(self, numbers):
542 """Set atomic numbers."""
543 self.set_array('numbers', numbers, int, ())
545 def get_atomic_numbers(self):
546 """Get integer array of atomic numbers."""
547 return self.arrays['numbers'].copy()
549 def get_chemical_symbols(self):
550 """Get list of chemical symbol strings.
552 Equivalent to ``list(atoms.symbols)``."""
553 return list(self.symbols)
555 def set_chemical_symbols(self, symbols):
556 """Set chemical symbols."""
557 self.set_array('numbers', symbols2numbers(symbols), int, ())
559 def get_chemical_formula(self, mode='hill', empirical=False):
560 """Get the chemical formula as a string based on the chemical symbols.
562 Parameters:
564 mode: str
565 There are four different modes available:
567 'all': The list of chemical symbols are contracted to a string,
568 e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'.
570 'reduce': The same as 'all' where repeated elements are contracted
571 to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to
572 'CH3OCH3'.
574 'hill': The list of chemical symbols are contracted to a string
575 following the Hill notation (alphabetical order with C and H
576 first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to
577 'H2O4S'. This is default.
579 'metal': The list of chemical symbols (alphabetical metals,
580 and alphabetical non-metals)
582 empirical, bool (optional, default=False)
583 Divide the symbol counts by their greatest common divisor to yield
584 an empirical formula. Only for mode `metal` and `hill`.
585 """
586 return self.symbols.get_chemical_formula(mode, empirical)
588 def set_tags(self, tags):
589 """Set tags for all atoms. If only one tag is supplied, it is
590 applied to all atoms."""
591 if isinstance(tags, int):
592 tags = [tags] * len(self)
593 self.set_array('tags', tags, int, ())
595 def get_tags(self):
596 """Get integer array of tags."""
597 if 'tags' in self.arrays:
598 return self.arrays['tags'].copy()
599 else:
600 return np.zeros(len(self), int)
602 def set_momenta(self, momenta, apply_constraint=True):
603 """Set momenta."""
604 if (apply_constraint and len(self.constraints) > 0 and
605 momenta is not None):
606 momenta = np.array(momenta) # modify a copy
607 for constraint in self.constraints:
608 if hasattr(constraint, 'adjust_momenta'):
609 constraint.adjust_momenta(self, momenta)
610 self.set_array('momenta', momenta, float, (3,))
612 def set_velocities(self, velocities):
613 """Set the momenta by specifying the velocities."""
614 self.set_momenta(self.get_masses()[:, np.newaxis] * velocities)
616 def get_momenta(self):
617 """Get array of momenta."""
618 if 'momenta' in self.arrays:
619 return self.arrays['momenta'].copy()
620 else:
621 return np.zeros((len(self), 3))
623 def set_masses(self, masses='defaults'):
624 """Set atomic masses in atomic mass units.
626 The array masses should contain a list of masses. In case
627 the masses argument is not given or for those elements of the
628 masses list that are None, standard values are set."""
630 if isinstance(masses, str):
631 if masses == 'defaults':
632 masses = atomic_masses[self.arrays['numbers']]
633 elif masses == 'most_common':
634 masses = atomic_masses_common[self.arrays['numbers']]
635 elif masses is None:
636 pass
637 elif not isinstance(masses, np.ndarray):
638 masses = list(masses)
639 for i, mass in enumerate(masses):
640 if mass is None:
641 masses[i] = atomic_masses[self.numbers[i]]
642 self.set_array('masses', masses, float, ())
644 def get_masses(self):
645 """Get array of masses in atomic mass units."""
646 if 'masses' in self.arrays:
647 return self.arrays['masses'].copy()
648 else:
649 return atomic_masses[self.arrays['numbers']]
651 def set_initial_magnetic_moments(self, magmoms=None):
652 """Set the initial magnetic moments.
654 Use either one or three numbers for every atom (collinear
655 or non-collinear spins)."""
657 if magmoms is None:
658 self.set_array('initial_magmoms', None)
659 else:
660 magmoms = np.asarray(magmoms)
661 self.set_array('initial_magmoms', magmoms, float,
662 magmoms.shape[1:])
664 def get_initial_magnetic_moments(self):
665 """Get array of initial magnetic moments."""
666 if 'initial_magmoms' in self.arrays:
667 return self.arrays['initial_magmoms'].copy()
668 else:
669 return np.zeros(len(self))
671 def get_magnetic_moments(self):
672 """Get calculated local magnetic moments."""
673 if self._calc is None:
674 raise RuntimeError('Atoms object has no calculator.')
675 return self._calc.get_magnetic_moments(self)
677 def get_magnetic_moment(self):
678 """Get calculated total magnetic moment."""
679 if self._calc is None:
680 raise RuntimeError('Atoms object has no calculator.')
681 return self._calc.get_magnetic_moment(self)
683 def set_initial_charges(self, charges=None):
684 """Set the initial charges."""
686 if charges is None:
687 self.set_array('initial_charges', None)
688 else:
689 self.set_array('initial_charges', charges, float, ())
691 def get_initial_charges(self):
692 """Get array of initial charges."""
693 if 'initial_charges' in self.arrays:
694 return self.arrays['initial_charges'].copy()
695 else:
696 return np.zeros(len(self))
698 def get_charges(self):
699 """Get calculated charges."""
700 if self._calc is None:
701 raise RuntimeError('Atoms object has no calculator.')
702 try:
703 return self._calc.get_charges(self)
704 except AttributeError:
705 from ase.calculators.calculator import PropertyNotImplementedError
706 raise PropertyNotImplementedError
708 def set_positions(self, newpositions, apply_constraint=True):
709 """Set positions, honoring any constraints. To ignore constraints,
710 use *apply_constraint=False*."""
711 if self.constraints and apply_constraint:
712 newpositions = np.array(newpositions, float)
713 for constraint in self.constraints:
714 constraint.adjust_positions(self, newpositions)
716 self.set_array('positions', newpositions, shape=(3,))
718 def get_positions(self, wrap=False, **wrap_kw):
719 """Get array of positions.
721 Parameters:
723 wrap: bool
724 wrap atoms back to the cell before returning positions
725 wrap_kw: (keyword=value) pairs
726 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
727 see :func:`ase.geometry.wrap_positions`
728 """
729 from ase.geometry import wrap_positions
730 if wrap:
731 if 'pbc' not in wrap_kw:
732 wrap_kw['pbc'] = self.pbc
733 return wrap_positions(self.positions, self.cell, **wrap_kw)
734 else:
735 return self.arrays['positions'].copy()
737 def get_potential_energy(self, force_consistent=False,
738 apply_constraint=True):
739 """Calculate potential energy.
741 Ask the attached calculator to calculate the potential energy and
742 apply constraints. Use *apply_constraint=False* to get the raw
743 forces.
745 When supported by the calculator, either the energy extrapolated
746 to zero Kelvin or the energy consistent with the forces (the free
747 energy) can be returned.
748 """
749 if self._calc is None:
750 raise RuntimeError('Atoms object has no calculator.')
751 if force_consistent:
752 energy = self._calc.get_potential_energy(
753 self, force_consistent=force_consistent)
754 else:
755 energy = self._calc.get_potential_energy(self)
756 if apply_constraint:
757 for constraint in self.constraints:
758 if hasattr(constraint, 'adjust_potential_energy'):
759 energy += constraint.adjust_potential_energy(self)
760 return energy
762 def get_properties(self, properties):
763 """This method is experimental; currently for internal use."""
764 # XXX Something about constraints.
765 if self._calc is None:
766 raise RuntimeError('Atoms object has no calculator.')
767 return self._calc.calculate_properties(self, properties)
769 def get_potential_energies(self):
770 """Calculate the potential energies of all the atoms.
772 Only available with calculators supporting per-atom energies
773 (e.g. classical potentials).
774 """
775 if self._calc is None:
776 raise RuntimeError('Atoms object has no calculator.')
777 return self._calc.get_potential_energies(self)
779 def get_kinetic_energy(self):
780 """Get the kinetic energy."""
781 momenta = self.arrays.get('momenta')
782 if momenta is None:
783 return 0.0
784 return 0.5 * np.vdot(momenta, self.get_velocities())
786 def get_velocities(self):
787 """Get array of velocities."""
788 momenta = self.get_momenta()
789 masses = self.get_masses()
790 return momenta / masses[:, np.newaxis]
792 def get_total_energy(self):
793 """Get the total energy - potential plus kinetic energy."""
794 return self.get_potential_energy() + self.get_kinetic_energy()
796 def get_forces(self, apply_constraint=True, md=False):
797 """Calculate atomic forces.
799 Ask the attached calculator to calculate the forces and apply
800 constraints. Use *apply_constraint=False* to get the raw
801 forces.
803 For molecular dynamics (md=True) we don't apply the constraint
804 to the forces but to the momenta. When holonomic constraints for
805 rigid linear triatomic molecules are present, ask the constraints
806 to redistribute the forces within each triple defined in the
807 constraints (required for molecular dynamics with this type of
808 constraints)."""
810 if self._calc is None:
811 raise RuntimeError('Atoms object has no calculator.')
812 forces = self._calc.get_forces(self)
814 if apply_constraint:
815 # We need a special md flag here because for MD we want
816 # to skip real constraints but include special "constraints"
817 # Like Hookean.
818 for constraint in self.constraints:
819 if md and hasattr(constraint, 'redistribute_forces_md'):
820 constraint.redistribute_forces_md(self, forces)
821 if not md or hasattr(constraint, 'adjust_potential_energy'):
822 constraint.adjust_forces(self, forces)
823 return forces
825 # Informs calculators (e.g. Asap) that ideal gas contribution is added here.
826 _ase_handles_dynamic_stress = True
828 def get_stress(self, voigt=True, apply_constraint=True,
829 include_ideal_gas=False):
830 """Calculate stress tensor.
832 Returns an array of the six independent components of the
833 symmetric stress tensor, in the traditional Voigt order
834 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt
835 order.
837 The ideal gas contribution to the stresses is added if the
838 atoms have momenta and ``include_ideal_gas`` is set to True.
839 """
841 if self._calc is None:
842 raise RuntimeError('Atoms object has no calculator.')
844 stress = self._calc.get_stress(self)
845 shape = stress.shape
847 if shape == (3, 3):
848 # Convert to the Voigt form before possibly applying
849 # constraints and adding the dynamic part of the stress
850 # (the "ideal gas contribution").
851 stress = full_3x3_to_voigt_6_stress(stress)
852 else:
853 assert shape == (6,)
855 if apply_constraint:
856 for constraint in self.constraints:
857 if hasattr(constraint, 'adjust_stress'):
858 constraint.adjust_stress(self, stress)
860 # Add ideal gas contribution, if applicable
861 if include_ideal_gas and self.has('momenta'):
862 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
863 p = self.get_momenta()
864 masses = self.get_masses()
865 invmass = 1.0 / masses
866 invvol = 1.0 / self.get_volume()
867 for alpha in range(3):
868 for beta in range(alpha, 3):
869 stress[stresscomp[alpha, beta]] -= (
870 p[:, alpha] * p[:, beta] * invmass).sum() * invvol
872 if voigt:
873 return stress
874 else:
875 return voigt_6_to_full_3x3_stress(stress)
877 def get_stresses(self, include_ideal_gas=False, voigt=True):
878 """Calculate the stress-tensor of all the atoms.
880 Only available with calculators supporting per-atom energies and
881 stresses (e.g. classical potentials). Even for such calculators
882 there is a certain arbitrariness in defining per-atom stresses.
884 The ideal gas contribution to the stresses is added if the
885 atoms have momenta and ``include_ideal_gas`` is set to True.
886 """
887 if self._calc is None:
888 raise RuntimeError('Atoms object has no calculator.')
889 stresses = self._calc.get_stresses(self)
891 # make sure `stresses` are in voigt form
892 if np.shape(stresses)[1:] == (3, 3):
893 stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses]
894 stresses = np.array(stresses_voigt)
896 # REMARK: The ideal gas contribution is intensive, i.e., the volume
897 # is divided out. We currently don't check if `stresses` are intensive
898 # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`.
899 # It might be good to check this here, but adds computational overhead.
901 if include_ideal_gas and self.has('momenta'):
902 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
903 if hasattr(self._calc, 'get_atomic_volumes'):
904 invvol = 1.0 / self._calc.get_atomic_volumes()
905 else:
906 invvol = self.get_global_number_of_atoms() / self.get_volume()
907 p = self.get_momenta()
908 invmass = 1.0 / self.get_masses()
909 for alpha in range(3):
910 for beta in range(alpha, 3):
911 stresses[:, stresscomp[alpha, beta]] -= (
912 p[:, alpha] * p[:, beta] * invmass * invvol)
913 if voigt:
914 return stresses
915 else:
916 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
917 return np.array(stresses_3x3)
919 def get_dipole_moment(self):
920 """Calculate the electric dipole moment for the atoms object.
922 Only available for calculators which has a get_dipole_moment()
923 method."""
925 if self._calc is None:
926 raise RuntimeError('Atoms object has no calculator.')
927 return self._calc.get_dipole_moment(self)
929 def copy(self):
930 """Return a copy."""
931 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
932 celldisp=self._celldisp.copy())
934 atoms.arrays = {}
935 for name, a in self.arrays.items():
936 atoms.arrays[name] = a.copy()
937 atoms.constraints = copy.deepcopy(self.constraints)
938 return atoms
940 def todict(self):
941 """For basic JSON (non-database) support."""
942 d = dict(self.arrays)
943 d['cell'] = np.asarray(self.cell)
944 d['pbc'] = self.pbc
945 if self._celldisp.any():
946 d['celldisp'] = self._celldisp
947 if self.constraints:
948 d['constraints'] = self.constraints
949 if self.info:
950 d['info'] = self.info
951 # Calculator... trouble.
952 return d
954 @classmethod
955 def fromdict(cls, dct):
956 """Rebuild atoms object from dictionary representation (todict)."""
957 dct = dct.copy()
958 kw = {name: dct.pop(name)
959 for name in ['numbers', 'positions', 'cell', 'pbc']}
960 constraints = dct.pop('constraints', None)
961 if constraints:
962 from ase.constraints import dict2constraint
963 constraints = [dict2constraint(d) for d in constraints]
965 info = dct.pop('info', None)
967 atoms = cls(constraint=constraints,
968 celldisp=dct.pop('celldisp', None),
969 info=info, **kw)
970 natoms = len(atoms)
972 # Some arrays are named differently from the atoms __init__ keywords.
973 # Also, there may be custom arrays. Hence we set them directly:
974 for name, arr in dct.items():
975 assert len(arr) == natoms, name
976 assert isinstance(arr, np.ndarray)
977 atoms.arrays[name] = arr
978 return atoms
980 def __len__(self):
981 return len(self.arrays['positions'])
983 @deprecated(
984 "Please use len(self) or, if your atoms are distributed, "
985 "self.get_global_number_of_atoms.",
986 category=FutureWarning,
987 )
988 def get_number_of_atoms(self):
989 """
990 .. deprecated:: 3.18.1
991 You probably want ``len(atoms)``. Or if your atoms are distributed,
992 use (and see) :func:`get_global_number_of_atoms()`.
993 """
994 return len(self)
996 def get_global_number_of_atoms(self):
997 """Returns the global number of atoms in a distributed-atoms parallel
998 simulation.
1000 DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING!
1002 Equivalent to len(atoms) in the standard ASE Atoms class. You should
1003 normally use len(atoms) instead. This function's only purpose is to
1004 make compatibility between ASE and Asap easier to maintain by having a
1005 few places in ASE use this function instead. It is typically only
1006 when counting the global number of degrees of freedom or in similar
1007 situations.
1008 """
1009 return len(self)
1011 def __repr__(self):
1012 tokens = []
1014 N = len(self)
1015 if N <= 60:
1016 symbols = self.get_chemical_formula('reduce')
1017 else:
1018 symbols = self.get_chemical_formula('hill')
1019 tokens.append(f"symbols='{symbols}'")
1021 if self.pbc.any() and not self.pbc.all():
1022 tokens.append(f'pbc={self.pbc.tolist()}')
1023 else:
1024 tokens.append(f'pbc={self.pbc[0]}')
1026 cell = self.cell
1027 if cell:
1028 if cell.orthorhombic:
1029 cell = cell.lengths().tolist()
1030 else:
1031 cell = cell.tolist()
1032 tokens.append(f'cell={cell}')
1034 for name in sorted(self.arrays):
1035 if name in ['numbers', 'positions']:
1036 continue
1037 tokens.append(f'{name}=...')
1039 if self.constraints:
1040 if len(self.constraints) == 1:
1041 constraint = self.constraints[0]
1042 else:
1043 constraint = self.constraints
1044 tokens.append(f'constraint={constraint!r}')
1046 if self._calc is not None:
1047 tokens.append('calculator={}(...)'
1048 .format(self._calc.__class__.__name__))
1050 return '{}({})'.format(self.__class__.__name__, ', '.join(tokens))
1052 def __add__(self, other):
1053 atoms = self.copy()
1054 atoms += other
1055 return atoms
1057 def extend(self, other):
1058 """Extend atoms object by appending atoms from *other*."""
1059 if isinstance(other, Atom):
1060 other = self.__class__([other])
1062 n1 = len(self)
1063 n2 = len(other)
1065 for name, a1 in self.arrays.items():
1066 a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype)
1067 a[:n1] = a1
1068 if name == 'masses':
1069 a2 = other.get_masses()
1070 else:
1071 a2 = other.arrays.get(name)
1072 if a2 is not None:
1073 a[n1:] = a2
1074 self.arrays[name] = a
1076 for name, a2 in other.arrays.items():
1077 if name in self.arrays:
1078 continue
1079 a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype)
1080 a[n1:] = a2
1081 if name == 'masses':
1082 a[:n1] = self.get_masses()[:n1]
1083 else:
1084 a[:n1] = 0
1086 self.set_array(name, a)
1088 def __iadd__(self, other):
1089 self.extend(other)
1090 return self
1092 def append(self, atom):
1093 """Append atom to end."""
1094 self.extend(self.__class__([atom]))
1096 def __iter__(self):
1097 for i in range(len(self)):
1098 yield self[i]
1100 def __getitem__(self, i):
1101 """Return a subset of the atoms.
1103 i -- scalar integer, list of integers, or slice object
1104 describing which atoms to return.
1106 If i is a scalar, return an Atom object. If i is a list or a
1107 slice, return an Atoms object with the same cell, pbc, and
1108 other associated info as the original Atoms object. The
1109 indices of the constraints will be shuffled so that they match
1110 the indexing in the subset returned.
1112 """
1114 if isinstance(i, numbers.Integral):
1115 natoms = len(self)
1116 if i < -natoms or i >= natoms:
1117 raise IndexError('Index out of range.')
1119 return Atom(atoms=self, index=i)
1120 elif not isinstance(i, slice):
1121 i = np.array(i)
1122 if len(i) == 0:
1123 i = np.array([], dtype=int)
1124 # if i is a mask
1125 if i.dtype == bool:
1126 if len(i) != len(self):
1127 raise IndexError('Length of mask {} must equal '
1128 'number of atoms {}'
1129 .format(len(i), len(self)))
1130 i = np.arange(len(self))[i]
1132 import copy
1134 conadd = []
1135 # Constraints need to be deepcopied, but only the relevant ones.
1136 for con in copy.deepcopy(self.constraints):
1137 try:
1138 con.index_shuffle(self, i)
1139 except (IndexError, NotImplementedError):
1140 pass
1141 else:
1142 conadd.append(con)
1144 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1145 # should be communicated to the slice as well
1146 celldisp=self._celldisp)
1147 # TODO: Do we need to shuffle indices in adsorbate_info too?
1149 atoms.arrays = {}
1150 for name, a in self.arrays.items():
1151 atoms.arrays[name] = a[i].copy()
1153 atoms.constraints = conadd
1154 return atoms
1156 def __delitem__(self, i):
1157 from ase.constraints import FixAtoms
1158 for c in self._constraints:
1159 if not isinstance(c, FixAtoms):
1160 raise RuntimeError('Remove constraint using set_constraint() '
1161 'before deleting atoms.')
1163 if isinstance(i, list) and len(i) > 0:
1164 # Make sure a list of booleans will work correctly and not be
1165 # interpreted at 0 and 1 indices.
1166 i = np.array(i)
1168 if len(self._constraints) > 0:
1169 n = len(self)
1170 i = np.arange(n)[i]
1171 if isinstance(i, int):
1172 i = [i]
1173 constraints = []
1174 for c in self._constraints:
1175 c = c.delete_atoms(i, n)
1176 if c is not None:
1177 constraints.append(c)
1178 self.constraints = constraints
1180 mask = np.ones(len(self), bool)
1181 mask[i] = False
1182 for name, a in self.arrays.items():
1183 self.arrays[name] = a[mask]
1185 def pop(self, i=-1):
1186 """Remove and return atom at index *i* (default last)."""
1187 atom = self[i]
1188 atom.cut_reference_to_atoms()
1189 del self[i]
1190 return atom
1192 def __imul__(self, m):
1193 """In-place repeat of atoms."""
1194 if isinstance(m, int):
1195 m = (m, m, m)
1197 for x, vec in zip(m, self.cell):
1198 if x != 1 and not vec.any():
1199 raise ValueError('Cannot repeat along undefined lattice '
1200 'vector')
1202 M = np.prod(m)
1203 n = len(self)
1205 for name, a in self.arrays.items():
1206 self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1))
1208 positions = self.arrays['positions']
1209 i0 = 0
1210 for m0 in range(m[0]):
1211 for m1 in range(m[1]):
1212 for m2 in range(m[2]):
1213 i1 = i0 + n
1214 positions[i0:i1] += np.dot((m0, m1, m2), self.cell)
1215 i0 = i1
1217 if self.constraints is not None:
1218 self.constraints = [c.repeat(m, n) for c in self.constraints]
1220 self.cell = np.array([m[c] * self.cell[c] for c in range(3)])
1222 return self
1224 def repeat(self, rep):
1225 """Create new repeated atoms object.
1227 The *rep* argument should be a sequence of three positive
1228 integers like *(2,3,1)* or a single integer (*r*) equivalent
1229 to *(r,r,r)*."""
1231 atoms = self.copy()
1232 atoms *= rep
1233 return atoms
1235 def __mul__(self, rep):
1236 return self.repeat(rep)
1238 def translate(self, displacement):
1239 """Translate atomic positions.
1241 The displacement argument can be a float an xyz vector or an
1242 nx3 array (where n is the number of atoms)."""
1244 self.arrays['positions'] += np.array(displacement)
1246 def center(self, vacuum=None, axis=(0, 1, 2), about=None):
1247 """Center atoms in unit cell.
1249 Centers the atoms in the unit cell, so there is the same
1250 amount of vacuum on all sides.
1252 vacuum: float (default: None)
1253 If specified adjust the amount of vacuum when centering.
1254 If vacuum=10.0 there will thus be 10 Angstrom of vacuum
1255 on each side.
1256 axis: int or sequence of ints
1257 Axis or axes to act on. Default: Act on all axes.
1258 about: float or array (default: None)
1259 If specified, center the atoms about <about>.
1260 I.e., about=(0., 0., 0.) (or just "about=0.", interpreted
1261 identically), to center about the origin.
1262 """
1264 # Find the orientations of the faces of the unit cell
1265 cell = self.cell.complete()
1266 dirs = np.zeros_like(cell)
1268 lengths = cell.lengths()
1269 for i in range(3):
1270 dirs[i] = np.cross(cell[i - 1], cell[i - 2])
1271 dirs[i] /= np.linalg.norm(dirs[i])
1272 if dirs[i] @ cell[i] < 0.0:
1273 dirs[i] *= -1
1275 if isinstance(axis, int):
1276 axes = (axis,)
1277 else:
1278 axes = axis
1280 # Now, decide how much each basis vector should be made longer
1281 pos = self.positions
1282 longer = np.zeros(3)
1283 shift = np.zeros(3)
1284 for i in axes:
1285 if len(pos):
1286 scalarprod = pos @ dirs[i]
1287 p0 = scalarprod.min()
1288 p1 = scalarprod.max()
1289 else:
1290 p0 = 0
1291 p1 = 0
1292 height = cell[i] @ dirs[i]
1293 if vacuum is not None:
1294 lng = (p1 - p0 + 2 * vacuum) - height
1295 else:
1296 lng = 0.0 # Do not change unit cell size!
1297 top = lng + height - p1
1298 shf = 0.5 * (top - p0)
1299 cosphi = cell[i] @ dirs[i] / lengths[i]
1300 longer[i] = lng / cosphi
1301 shift[i] = shf / cosphi
1303 # Now, do it!
1304 translation = np.zeros(3)
1305 for i in axes:
1306 nowlen = lengths[i]
1307 if vacuum is not None:
1308 self.cell[i] = cell[i] * (1 + longer[i] / nowlen)
1309 translation += shift[i] * cell[i] / nowlen
1311 # We calculated translations using the completed cell,
1312 # so directions without cell vectors will have been centered
1313 # along a "fake" vector of length 1.
1314 # Therefore, we adjust by -0.5:
1315 if not any(self.cell[i]):
1316 translation[i] -= 0.5
1318 # Optionally, translate to center about a point in space.
1319 if about is not None:
1320 for vector in self.cell:
1321 translation -= vector / 2.0
1322 translation += about
1324 self.positions += translation
1326 def get_center_of_mass(self, scaled=False, indices=None):
1327 """Get the center of mass.
1329 Parameters
1330 ----------
1331 scaled : bool
1332 If True, the center of mass in scaled coordinates is returned.
1333 indices : list | slice | str, default: None
1334 If specified, the center of mass of a subset of atoms is returned.
1335 """
1336 if indices is None:
1337 indices = slice(None)
1338 elif isinstance(indices, str):
1339 indices = string2index(indices)
1341 masses = self.get_masses()[indices]
1342 com = masses @ self.positions[indices] / masses.sum()
1343 if scaled:
1344 return self.cell.scaled_positions(com)
1345 return com # Cartesian coordinates
1347 def set_center_of_mass(self, com, scaled=False):
1348 """Set the center of mass.
1350 If scaled=True the center of mass is expected in scaled coordinates.
1351 Constraints are considered for scaled=False.
1352 """
1353 old_com = self.get_center_of_mass(scaled=scaled)
1354 difference = com - old_com
1355 if scaled:
1356 self.set_scaled_positions(self.get_scaled_positions() + difference)
1357 else:
1358 self.set_positions(self.get_positions() + difference)
1360 def get_moments_of_inertia(self, vectors=False):
1361 """Get the moments of inertia along the principal axes.
1363 The three principal moments of inertia are computed from the
1364 eigenvalues of the symmetric inertial tensor. Periodic boundary
1365 conditions are ignored. Units of the moments of inertia are
1366 amu*angstrom**2.
1367 """
1368 com = self.get_center_of_mass()
1369 positions = self.get_positions()
1370 positions -= com # translate center of mass to origin
1371 masses = self.get_masses()
1373 # Initialize elements of the inertial tensor
1374 I11 = I22 = I33 = I12 = I13 = I23 = 0.0
1375 for i in range(len(self)):
1376 x, y, z = positions[i]
1377 m = masses[i]
1379 I11 += m * (y ** 2 + z ** 2)
1380 I22 += m * (x ** 2 + z ** 2)
1381 I33 += m * (x ** 2 + y ** 2)
1382 I12 += -m * x * y
1383 I13 += -m * x * z
1384 I23 += -m * y * z
1386 Itensor = np.array([[I11, I12, I13],
1387 [I12, I22, I23],
1388 [I13, I23, I33]])
1390 evals, evecs = np.linalg.eigh(Itensor)
1391 if vectors:
1392 return evals, evecs.transpose()
1393 else:
1394 return evals
1396 def get_angular_momentum(self):
1397 """Get total angular momentum with respect to the center of mass."""
1398 com = self.get_center_of_mass()
1399 positions = self.get_positions()
1400 positions -= com # translate center of mass to origin
1401 return np.cross(positions, self.get_momenta()).sum(0)
1403 def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False):
1404 """Rotate atoms based on a vector and an angle, or two vectors.
1406 Parameters:
1408 a = None:
1409 Angle that the atoms is rotated around the vector 'v'. 'a'
1410 can also be a vector and then 'a' is rotated
1411 into 'v'.
1413 v:
1414 Vector to rotate the atoms around. Vectors can be given as
1415 strings: 'x', '-x', 'y', ... .
1417 center = (0, 0, 0):
1418 The center is kept fixed under the rotation. Use 'COM' to fix
1419 the center of mass, 'COP' to fix the center of positions or
1420 'COU' to fix the center of cell.
1422 rotate_cell = False:
1423 If true the cell is also rotated.
1425 Examples:
1427 Rotate 90 degrees around the z-axis, so that the x-axis is
1428 rotated into the y-axis:
1430 >>> atoms = Atoms()
1431 >>> atoms.rotate(90, 'z')
1432 >>> atoms.rotate(90, (0, 0, 1))
1433 >>> atoms.rotate(-90, '-z')
1434 >>> atoms.rotate('x', 'y')
1435 >>> atoms.rotate((1, 0, 0), (0, 1, 0))
1436 """
1438 if not isinstance(a, numbers.Real):
1439 a, v = v, a
1441 norm = np.linalg.norm
1442 v = string2vector(v)
1444 normv = norm(v)
1446 if normv == 0.0:
1447 raise ZeroDivisionError('Cannot rotate: norm(v) == 0')
1449 if isinstance(a, numbers.Real):
1450 a *= pi / 180
1451 v /= normv
1452 c = cos(a)
1453 s = sin(a)
1454 else:
1455 v2 = string2vector(a)
1456 v /= normv
1457 normv2 = np.linalg.norm(v2)
1458 if normv2 == 0:
1459 raise ZeroDivisionError('Cannot rotate: norm(a) == 0')
1460 v2 /= norm(v2)
1461 c = np.dot(v, v2)
1462 v = np.cross(v, v2)
1463 s = norm(v)
1464 # In case *v* and *a* are parallel, np.cross(v, v2) vanish
1465 # and can't be used as a rotation axis. However, in this
1466 # case any rotation axis perpendicular to v2 will do.
1467 eps = 1e-7
1468 if s < eps:
1469 v = np.cross((0, 0, 1), v2)
1470 if norm(v) < eps:
1471 v = np.cross((1, 0, 0), v2)
1472 assert norm(v) >= eps
1473 elif s > 0:
1474 v /= s
1476 center = self._centering_as_array(center)
1478 p = self.arrays['positions'] - center
1479 self.arrays['positions'][:] = (c * p -
1480 np.cross(p, s * v) +
1481 np.outer(np.dot(p, v), (1.0 - c) * v) +
1482 center)
1483 if rotate_cell:
1484 rotcell = self.get_cell()
1485 rotcell[:] = (c * rotcell -
1486 np.cross(rotcell, s * v) +
1487 np.outer(np.dot(rotcell, v), (1.0 - c) * v))
1488 self.set_cell(rotcell)
1490 def _centering_as_array(self, center):
1491 if isinstance(center, str):
1492 if center.lower() == 'com':
1493 center = self.get_center_of_mass()
1494 elif center.lower() == 'cop':
1495 center = self.get_positions().mean(axis=0)
1496 elif center.lower() == 'cou':
1497 center = self.get_cell().sum(axis=0) / 2
1498 else:
1499 raise ValueError('Cannot interpret center')
1500 else:
1501 center = np.array(center, float)
1502 return center
1504 def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)):
1505 """Rotate atoms via Euler angles (in degrees).
1507 See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation.
1509 Parameters:
1511 center :
1512 The point to rotate about. A sequence of length 3 with the
1513 coordinates, or 'COM' to select the center of mass, 'COP' to
1514 select center of positions or 'COU' to select center of cell.
1515 phi :
1516 The 1st rotation angle around the z axis.
1517 theta :
1518 Rotation around the x axis.
1519 psi :
1520 2nd rotation around the z axis.
1522 """
1523 center = self._centering_as_array(center)
1525 phi *= pi / 180
1526 theta *= pi / 180
1527 psi *= pi / 180
1529 # First move the molecule to the origin In contrast to MATLAB,
1530 # numpy broadcasts the smaller array to the larger row-wise,
1531 # so there is no need to play with the Kronecker product.
1532 rcoords = self.positions - center
1533 # First Euler rotation about z in matrix form
1534 D = np.array(((cos(phi), sin(phi), 0.),
1535 (-sin(phi), cos(phi), 0.),
1536 (0., 0., 1.)))
1537 # Second Euler rotation about x:
1538 C = np.array(((1., 0., 0.),
1539 (0., cos(theta), sin(theta)),
1540 (0., -sin(theta), cos(theta))))
1541 # Third Euler rotation, 2nd rotation about z:
1542 B = np.array(((cos(psi), sin(psi), 0.),
1543 (-sin(psi), cos(psi), 0.),
1544 (0., 0., 1.)))
1545 # Total Euler rotation
1546 A = np.dot(B, np.dot(C, D))
1547 # Do the rotation
1548 rcoords = np.dot(A, np.transpose(rcoords))
1549 # Move back to the rotation point
1550 self.positions = np.transpose(rcoords) + center
1552 def get_dihedral(self, a0, a1, a2, a3, mic=False):
1553 """Calculate dihedral angle.
1555 Calculate dihedral angle (in degrees) between the vectors a0->a1
1556 and a2->a3.
1558 Use mic=True to use the Minimum Image Convention and calculate the
1559 angle across periodic boundaries.
1560 """
1561 return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0]
1563 def get_dihedrals(self, indices, mic=False):
1564 """Calculate dihedral angles.
1566 Calculate dihedral angles (in degrees) between the list of vectors
1567 a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices.
1569 Use mic=True to use the Minimum Image Convention and calculate the
1570 angles across periodic boundaries.
1571 """
1572 from ase.geometry import get_dihedrals
1574 indices = np.array(indices)
1575 assert indices.shape[1] == 4
1577 a0s = self.positions[indices[:, 0]]
1578 a1s = self.positions[indices[:, 1]]
1579 a2s = self.positions[indices[:, 2]]
1580 a3s = self.positions[indices[:, 3]]
1582 # vectors 0->1, 1->2, 2->3
1583 v0 = a1s - a0s
1584 v1 = a2s - a1s
1585 v2 = a3s - a2s
1587 cell = None
1588 pbc = None
1590 if mic:
1591 cell = self.cell
1592 pbc = self.pbc
1594 return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc)
1596 def _masked_rotate(self, center, axis, diff, mask):
1597 # do rotation of subgroup by copying it to temporary atoms object
1598 # and then rotating that
1599 #
1600 # recursive object definition might not be the most elegant thing,
1601 # more generally useful might be a rotation function with a mask?
1602 group = self.__class__()
1603 for i in range(len(self)):
1604 if mask[i]:
1605 group += self[i]
1606 group.translate(-center)
1607 group.rotate(diff * 180 / pi, axis)
1608 group.translate(center)
1609 # set positions in original atoms object
1610 j = 0
1611 for i in range(len(self)):
1612 if mask[i]:
1613 self.positions[i] = group[j].position
1614 j += 1
1616 def set_dihedral(self, a1, a2, a3, a4, angle,
1617 mask=None, indices=None):
1618 """Set the dihedral angle (degrees) between vectors a1->a2 and
1619 a3->a4 by changing the atom indexed by a4.
1621 If mask is not None, all the atoms described in mask
1622 (read: the entire subgroup) are moved. Alternatively to the mask,
1623 the indices of the atoms to be rotated can be supplied. If both
1624 *mask* and *indices* are given, *indices* overwrites *mask*.
1626 **Important**: If *mask* or *indices* is given and does not contain
1627 *a4*, *a4* will NOT be moved. In most cases you therefore want
1628 to include *a4* in *mask*/*indices*.
1630 Example: the following defines a very crude
1631 ethane-like molecule and twists one half of it by 30 degrees.
1633 >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0],
1634 ... [1, 0, 0], [2, 1, 0], [2, -1, 0]])
1635 >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1])
1636 """
1638 angle *= pi / 180
1640 # if not provided, set mask to the last atom in the
1641 # dihedral description
1642 if mask is None and indices is None:
1643 mask = np.zeros(len(self))
1644 mask[a4] = 1
1645 elif indices is not None:
1646 mask = [index in indices for index in range(len(self))]
1648 # compute necessary in dihedral change, from current value
1649 current = self.get_dihedral(a1, a2, a3, a4) * pi / 180
1650 diff = angle - current
1651 axis = self.positions[a3] - self.positions[a2]
1652 center = self.positions[a3]
1653 self._masked_rotate(center, axis, diff, mask)
1655 def rotate_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None):
1656 """Rotate dihedral angle.
1658 Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a
1659 predefined dihedral angle, starting from its current configuration.
1660 """
1661 start = self.get_dihedral(a1, a2, a3, a4)
1662 self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices)
1664 def get_angle(self, a1, a2, a3, mic=False):
1665 """Get angle formed by three atoms.
1667 Calculate angle in degrees between the vectors a2->a1 and
1668 a2->a3.
1670 Use mic=True to use the Minimum Image Convention and calculate the
1671 angle across periodic boundaries.
1672 """
1673 return self.get_angles([[a1, a2, a3]], mic=mic)[0]
1675 def get_angles(self, indices, mic=False):
1676 """Get angle formed by three atoms for multiple groupings.
1678 Calculate angle in degrees between vectors between atoms a2->a1
1679 and a2->a3, where a1, a2, and a3 are in each row of indices.
1681 Use mic=True to use the Minimum Image Convention and calculate
1682 the angle across periodic boundaries.
1683 """
1684 from ase.geometry import get_angles
1686 indices = np.array(indices)
1687 assert indices.shape[1] == 3
1689 a1s = self.positions[indices[:, 0]]
1690 a2s = self.positions[indices[:, 1]]
1691 a3s = self.positions[indices[:, 2]]
1693 v12 = a1s - a2s
1694 v32 = a3s - a2s
1696 cell = None
1697 pbc = None
1699 if mic:
1700 cell = self.cell
1701 pbc = self.pbc
1703 return get_angles(v12, v32, cell=cell, pbc=pbc)
1705 def set_angle(self, a1, a2=None, a3=None, angle=None, mask=None,
1706 indices=None, add=False):
1707 """Set angle (in degrees) formed by three atoms.
1709 Sets the angle between vectors *a2*->*a1* and *a2*->*a3*.
1711 If *add* is `True`, the angle will be changed by the value given.
1713 Same usage as in :meth:`ase.Atoms.set_dihedral`.
1714 If *mask* and *indices*
1715 are given, *indices* overwrites *mask*. If *mask* and *indices*
1716 are not set, only *a3* is moved."""
1718 if any(a is None for a in [a2, a3, angle]):
1719 raise ValueError('a2, a3, and angle must not be None')
1721 # If not provided, set mask to the last atom in the angle description
1722 if mask is None and indices is None:
1723 mask = np.zeros(len(self))
1724 mask[a3] = 1
1725 elif indices is not None:
1726 mask = [index in indices for index in range(len(self))]
1728 if add:
1729 diff = angle
1730 else:
1731 # Compute necessary in angle change, from current value
1732 diff = angle - self.get_angle(a1, a2, a3)
1734 diff *= pi / 180
1735 # Do rotation of subgroup by copying it to temporary atoms object and
1736 # then rotating that
1737 v10 = self.positions[a1] - self.positions[a2]
1738 v12 = self.positions[a3] - self.positions[a2]
1739 v10 /= np.linalg.norm(v10)
1740 v12 /= np.linalg.norm(v12)
1741 axis = np.cross(v10, v12)
1742 center = self.positions[a2]
1743 self._masked_rotate(center, axis, diff, mask)
1745 def rattle(self, stdev=0.001, seed=None, rng=None):
1746 """Randomly displace atoms.
1748 This method adds random displacements to the atomic positions,
1749 taking a possible constraint into account. The random numbers are
1750 drawn from a normal distribution of standard deviation stdev.
1752 By default, the random number generator always uses the same seed (42)
1753 for repeatability. You can provide your own seed (an integer), or if you
1754 want the randomness to be different each time you run a script, then
1755 provide `rng=numpy.random`. For a parallel calculation, it is important
1756 to use the same seed on all processors! """
1758 if seed is not None and rng is not None:
1759 raise ValueError('Please do not provide both seed and rng.')
1761 if rng is None:
1762 if seed is None:
1763 seed = 42
1764 rng = np.random.RandomState(seed)
1765 positions = self.arrays['positions']
1766 self.set_positions(positions +
1767 rng.normal(scale=stdev, size=positions.shape))
1769 def get_distance(self, a0, a1, mic=False, vector=False):
1770 """Return distance between two atoms.
1772 Use mic=True to use the Minimum Image Convention.
1773 vector=True gives the distance vector (from a0 to a1).
1774 """
1775 return self.get_distances(a0, [a1], mic=mic, vector=vector)[0]
1777 def get_distances(self, a, indices, mic=False, vector=False):
1778 """Return distances of atom No.i with a list of atoms.
1780 Use mic=True to use the Minimum Image Convention.
1781 vector=True gives the distance vector (from a to self[indices]).
1782 """
1783 from ase.geometry import get_distances
1785 R = self.arrays['positions']
1786 p1 = [R[a]]
1787 p2 = R[indices]
1789 cell = None
1790 pbc = None
1792 if mic:
1793 cell = self.cell
1794 pbc = self.pbc
1796 D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc)
1798 if vector:
1799 D.shape = (-1, 3)
1800 return D
1801 else:
1802 D_len.shape = (-1,)
1803 return D_len
1805 def get_all_distances(self, mic=False, vector=False):
1806 """Return distances of all of the atoms with all of the atoms.
1808 Use mic=True to use the Minimum Image Convention.
1809 """
1810 from ase.geometry import get_distances
1812 R = self.arrays['positions']
1814 cell = None
1815 pbc = None
1817 if mic:
1818 cell = self.cell
1819 pbc = self.pbc
1821 D, D_len = get_distances(R, cell=cell, pbc=pbc)
1823 if vector:
1824 return D
1825 else:
1826 return D_len
1828 def set_distance(self, a0, a1, distance, fix=0.5, mic=False,
1829 mask=None, indices=None, add=False, factor=False):
1830 """Set the distance between two atoms.
1832 Set the distance between atoms *a0* and *a1* to *distance*.
1833 By default, the center of the two atoms will be fixed. Use
1834 *fix=0* to fix the first atom, *fix=1* to fix the second
1835 atom and *fix=0.5* (default) to fix the center of the bond.
1837 If *mask* or *indices* are set (*mask* overwrites *indices*),
1838 only the atoms defined there are moved
1839 (see :meth:`ase.Atoms.set_dihedral`).
1841 When *add* is true, the distance is changed by the value given.
1842 In combination
1843 with *factor* True, the value given is a factor scaling the distance.
1845 It is assumed that the atoms in *mask*/*indices* move together
1846 with *a1*. If *fix=1*, only *a0* will therefore be moved."""
1847 from ase.geometry import find_mic
1849 if a0 % len(self) == a1 % len(self):
1850 raise ValueError('a0 and a1 must not be the same')
1852 if add:
1853 oldDist = self.get_distance(a0, a1, mic=mic)
1854 if factor:
1855 newDist = oldDist * distance
1856 else:
1857 newDist = oldDist + distance
1858 self.set_distance(a0, a1, newDist, fix=fix, mic=mic,
1859 mask=mask, indices=indices, add=False,
1860 factor=False)
1861 return
1863 R = self.arrays['positions']
1864 D = np.array([R[a1] - R[a0]])
1866 if mic:
1867 D, D_len = find_mic(D, self.cell, self.pbc)
1868 else:
1869 D_len = np.array([np.sqrt((D**2).sum())])
1870 x = 1.0 - distance / D_len[0]
1872 if mask is None and indices is None:
1873 indices = [a0, a1]
1874 elif mask:
1875 indices = [i for i in range(len(self)) if mask[i]]
1877 for i in indices:
1878 if i == a0:
1879 R[a0] += (x * fix) * D[0]
1880 else:
1881 R[i] -= (x * (1.0 - fix)) * D[0]
1883 def get_scaled_positions(self, wrap=True):
1884 """Get positions relative to unit cell.
1886 If wrap is True, atoms outside the unit cell will be wrapped into
1887 the cell in those directions with periodic boundary conditions
1888 so that the scaled coordinates are between zero and one.
1890 If any cell vectors are zero, the corresponding coordinates
1891 are evaluated as if the cell were completed using
1892 ``cell.complete()``. This means coordinates will be Cartesian
1893 as long as the non-zero cell vectors span a Cartesian axis or
1894 plane."""
1896 fractional = self.cell.scaled_positions(self.positions)
1898 if wrap:
1899 for i, periodic in enumerate(self.pbc):
1900 if periodic:
1901 # Yes, we need to do it twice.
1902 # See the scaled_positions.py test.
1903 fractional[:, i] %= 1.0
1904 fractional[:, i] %= 1.0
1906 return fractional
1908 def set_scaled_positions(self, scaled):
1909 """Set positions relative to unit cell."""
1910 self.positions[:] = self.cell.cartesian_positions(scaled)
1912 def wrap(self, **wrap_kw):
1913 """Wrap positions to unit cell.
1915 Parameters:
1917 wrap_kw: (keyword=value) pairs
1918 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
1919 see :func:`ase.geometry.wrap_positions`
1920 """
1922 if 'pbc' not in wrap_kw:
1923 wrap_kw['pbc'] = self.pbc
1925 self.positions[:] = self.get_positions(wrap=True, **wrap_kw)
1927 def get_temperature(self):
1928 """Get the temperature in Kelvin."""
1929 dof = len(self) * 3
1930 for constraint in self._constraints:
1931 dof -= constraint.get_removed_dof(self)
1932 ekin = self.get_kinetic_energy()
1933 return 2 * ekin / (dof * units.kB)
1935 def __eq__(self, other):
1936 """Check for identity of two atoms objects.
1938 Identity means: same positions, atomic numbers, unit cell and
1939 periodic boundary conditions."""
1940 if not isinstance(other, Atoms):
1941 return False
1942 a = self.arrays
1943 b = other.arrays
1944 return (len(self) == len(other) and
1945 (a['positions'] == b['positions']).all() and
1946 (a['numbers'] == b['numbers']).all() and
1947 (self.cell == other.cell).all() and
1948 (self.pbc == other.pbc).all())
1950 def __ne__(self, other):
1951 """Check if two atoms objects are not equal.
1953 Any differences in positions, atomic numbers, unit cell or
1954 periodic boundary condtions make atoms objects not equal.
1955 """
1956 eq = self.__eq__(other)
1957 if eq is NotImplemented:
1958 return eq
1959 else:
1960 return not eq
1962 # @deprecated('Please use atoms.cell.volume')
1963 # We kind of want to deprecate this, but the ValueError behaviour
1964 # might be desirable. Should we do this?
1965 def get_volume(self):
1966 """Get volume of unit cell."""
1967 if self.cell.rank != 3:
1968 raise ValueError(
1969 'You have {} lattice vectors: volume not defined'
1970 .format(self.cell.rank))
1971 return self.cell.volume
1973 def _get_positions(self):
1974 """Return reference to positions-array for in-place manipulations."""
1975 return self.arrays['positions']
1977 def _set_positions(self, pos):
1978 """Set positions directly, bypassing constraints."""
1979 self.arrays['positions'][:] = pos
1981 positions = property(_get_positions, _set_positions,
1982 doc='Attribute for direct ' +
1983 'manipulation of the positions.')
1985 def _get_atomic_numbers(self):
1986 """Return reference to atomic numbers for in-place
1987 manipulations."""
1988 return self.arrays['numbers']
1990 numbers = property(_get_atomic_numbers, set_atomic_numbers,
1991 doc='Attribute for direct ' +
1992 'manipulation of the atomic numbers.')
1994 @property
1995 def cell(self):
1996 """The :class:`ase.cell.Cell` for direct manipulation."""
1997 return self._cellobj
1999 @cell.setter
2000 def cell(self, cell):
2001 cell = Cell.ascell(cell)
2002 self._cellobj[:] = cell
2004 def write(self, filename, format=None, **kwargs):
2005 """Write atoms object to a file.
2007 see ase.io.write for formats.
2008 kwargs are passed to ase.io.write.
2009 """
2010 from ase.io import write
2011 write(filename, self, format, **kwargs)
2013 def iterimages(self):
2014 yield self
2016 def __ase_optimizable__(self):
2017 from ase.optimize.optimize import OptimizableAtoms
2018 return OptimizableAtoms(self)
2020 def edit(self):
2021 """Modify atoms interactively through ASE's GUI viewer.
2023 Conflicts leading to undesirable behaviour might arise
2024 when matplotlib has been pre-imported with certain
2025 incompatible backends and while trying to use the
2026 plot feature inside the interactive GUI. To circumvent,
2027 please set matplotlib.use('gtk') before calling this
2028 method.
2029 """
2030 from ase.gui.gui import GUI
2031 from ase.gui.images import Images
2032 images = Images([self])
2033 gui = GUI(images)
2034 gui.run()
2037def string2vector(v):
2038 if isinstance(v, str):
2039 if v[0] == '-':
2040 return -string2vector(v[1:])
2041 w = np.zeros(3)
2042 w['xyz'.index(v)] = 1.0
2043 return w
2044 return np.array(v, float)
2047def default(data, dflt):
2048 """Helper function for setting default values."""
2049 if data is None:
2050 return None
2051 elif isinstance(data, (list, tuple)):
2052 newdata = []
2053 allnone = True
2054 for x in data:
2055 if x is None:
2056 newdata.append(dflt)
2057 else:
2058 newdata.append(x)
2059 allnone = False
2060 if allnone:
2061 return None
2062 return newdata
2063 else:
2064 return data