Coverage for /builds/alexhroom/ase/ase/atoms.py: 93.89%
965 statements
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-05 14:37 +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 get_number_of_degrees_of_freedom(self):
349 """Calculate the number of degrees of freedom in the system."""
350 return len(self) * 3 - sum(
351 c.get_removed_dof(self) for c in self._constraints
352 )
354 def set_cell(self, cell, scale_atoms=False, apply_constraint=True):
355 """Set unit cell vectors.
357 Parameters:
359 cell: 3x3 matrix or length 3 or 6 vector
360 Unit cell. A 3x3 matrix (the three unit cell vectors) or
361 just three numbers for an orthorhombic cell. Another option is
362 6 numbers, which describes unit cell with lengths of unit cell
363 vectors and with angles between them (in degrees), in following
364 order: [len(a), len(b), len(c), angle(b,c), angle(a,c),
365 angle(a,b)]. First vector will lie in x-direction, second in
366 xy-plane, and the third one in z-positive subspace.
367 scale_atoms: bool
368 Fix atomic positions or move atoms with the unit cell?
369 Default behavior is to *not* move the atoms (scale_atoms=False).
370 apply_constraint: bool
371 Whether to apply constraints to the given cell.
373 Examples:
375 Two equivalent ways to define an orthorhombic cell:
377 >>> atoms = Atoms('He')
378 >>> a, b, c = 7, 7.5, 8
379 >>> atoms.set_cell([a, b, c])
380 >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)])
382 FCC unit cell:
384 >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)])
386 Hexagonal unit cell:
388 >>> atoms.set_cell([a, a, c, 90, 90, 120])
390 Rhombohedral unit cell:
392 >>> alpha = 77
393 >>> atoms.set_cell([a, a, a, alpha, alpha, alpha])
394 """
396 # Override pbcs if and only if given a Cell object:
397 cell = Cell.new(cell)
399 # XXX not working well during initialize due to missing _constraints
400 if apply_constraint and hasattr(self, '_constraints'):
401 for constraint in self.constraints:
402 if hasattr(constraint, 'adjust_cell'):
403 constraint.adjust_cell(self, cell)
405 if scale_atoms:
406 M = np.linalg.solve(self.cell.complete(), cell.complete())
407 self.positions[:] = np.dot(self.positions, M)
409 self.cell[:] = cell
411 def set_celldisp(self, celldisp):
412 """Set the unit cell displacement vectors."""
413 celldisp = np.array(celldisp, float)
414 self._celldisp = celldisp
416 def get_celldisp(self):
417 """Get the unit cell displacement vectors."""
418 return self._celldisp.copy()
420 def get_cell(self, complete=False):
421 """Get the three unit cell vectors as a `class`:ase.cell.Cell` object.
423 The Cell object resembles a 3x3 ndarray, and cell[i, j]
424 is the jth Cartesian coordinate of the ith cell vector."""
425 if complete:
426 cell = self.cell.complete()
427 else:
428 cell = self.cell.copy()
430 return cell
432 @deprecated('Please use atoms.cell.cellpar() instead', DeprecationWarning)
433 def get_cell_lengths_and_angles(self):
434 """Get unit cell parameters. Sequence of 6 numbers.
436 First three are unit cell vector lengths and second three
437 are angles between them::
439 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
441 in degrees.
443 .. deprecated:: 3.21.0
444 Please use ``atoms.cell.cellpar()`` instead
445 """
446 return self.cell.cellpar()
448 @deprecated('Please use atoms.cell.reciprocal()', DeprecationWarning)
449 def get_reciprocal_cell(self):
450 """Get the three reciprocal lattice vectors as a 3x3 ndarray.
452 Note that the commonly used factor of 2 pi for Fourier
453 transforms is not included here.
455 .. deprecated:: 3.21.0
456 Please use ``atoms.cell.reciprocal()``
457 """
458 return self.cell.reciprocal()
460 @property
461 def pbc(self):
462 """Reference to pbc-flags for in-place manipulations."""
463 return self._pbc
465 @pbc.setter
466 def pbc(self, pbc):
467 self._pbc[:] = pbc
469 def set_pbc(self, pbc):
470 """Set periodic boundary condition flags."""
471 self.pbc = pbc
473 def get_pbc(self):
474 """Get periodic boundary condition flags."""
475 return self.pbc.copy()
477 def new_array(self, name, a, dtype=None, shape=None):
478 """Add new array.
480 If *shape* is not *None*, the shape of *a* will be checked."""
482 if dtype is not None:
483 a = np.array(a, dtype, order='C')
484 if len(a) == 0 and shape is not None:
485 a.shape = (-1,) + shape
486 else:
487 if not a.flags['C_CONTIGUOUS']:
488 a = np.ascontiguousarray(a)
489 else:
490 a = a.copy()
492 if name in self.arrays:
493 raise RuntimeError(f'Array {name} already present')
495 for b in self.arrays.values():
496 if len(a) != len(b):
497 raise ValueError('Array "%s" has wrong length: %d != %d.' %
498 (name, len(a), len(b)))
499 break
501 if shape is not None and a.shape[1:] != shape:
502 raise ValueError(
503 f'Array "{name}" has wrong shape {a.shape} != '
504 f'{(a.shape[0:1] + shape)}.')
506 self.arrays[name] = a
508 def get_array(self, name, copy=True):
509 """Get an array.
511 Returns a copy unless the optional argument copy is false.
512 """
513 if copy:
514 return self.arrays[name].copy()
515 else:
516 return self.arrays[name]
518 def set_array(self, name, a, dtype=None, shape=None):
519 """Update array.
521 If *shape* is not *None*, the shape of *a* will be checked.
522 If *a* is *None*, then the array is deleted."""
524 b = self.arrays.get(name)
525 if b is None:
526 if a is not None:
527 self.new_array(name, a, dtype, shape)
528 else:
529 if a is None:
530 del self.arrays[name]
531 else:
532 a = np.asarray(a)
533 if a.shape != b.shape:
534 raise ValueError(
535 f'Array "{name}" has wrong shape '
536 f'{a.shape} != {b.shape}.')
537 b[:] = a
539 def has(self, name):
540 """Check for existence of array.
542 name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms',
543 'initial_charges'."""
544 # XXX extend has to calculator properties
545 return name in self.arrays
547 def set_atomic_numbers(self, numbers):
548 """Set atomic numbers."""
549 self.set_array('numbers', numbers, int, ())
551 def get_atomic_numbers(self):
552 """Get integer array of atomic numbers."""
553 return self.arrays['numbers'].copy()
555 def get_chemical_symbols(self):
556 """Get list of chemical symbol strings.
558 Equivalent to ``list(atoms.symbols)``."""
559 return list(self.symbols)
561 def set_chemical_symbols(self, symbols):
562 """Set chemical symbols."""
563 self.set_array('numbers', symbols2numbers(symbols), int, ())
565 def get_chemical_formula(self, mode='hill', empirical=False):
566 """Get the chemical formula as a string based on the chemical symbols.
568 Parameters:
570 mode: str
571 There are four different modes available:
573 'all': The list of chemical symbols are contracted to a string,
574 e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'.
576 'reduce': The same as 'all' where repeated elements are contracted
577 to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to
578 'CH3OCH3'.
580 'hill': The list of chemical symbols are contracted to a string
581 following the Hill notation (alphabetical order with C and H
582 first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to
583 'H2O4S'. This is default.
585 'metal': The list of chemical symbols (alphabetical metals,
586 and alphabetical non-metals)
588 empirical, bool (optional, default=False)
589 Divide the symbol counts by their greatest common divisor to yield
590 an empirical formula. Only for mode `metal` and `hill`.
591 """
592 return self.symbols.get_chemical_formula(mode, empirical)
594 def set_tags(self, tags):
595 """Set tags for all atoms. If only one tag is supplied, it is
596 applied to all atoms."""
597 if isinstance(tags, int):
598 tags = [tags] * len(self)
599 self.set_array('tags', tags, int, ())
601 def get_tags(self):
602 """Get integer array of tags."""
603 if 'tags' in self.arrays:
604 return self.arrays['tags'].copy()
605 else:
606 return np.zeros(len(self), int)
608 def set_momenta(self, momenta, apply_constraint=True):
609 """Set momenta."""
610 if (apply_constraint and len(self.constraints) > 0 and
611 momenta is not None):
612 momenta = np.array(momenta) # modify a copy
613 for constraint in self.constraints:
614 if hasattr(constraint, 'adjust_momenta'):
615 constraint.adjust_momenta(self, momenta)
616 self.set_array('momenta', momenta, float, (3,))
618 def set_velocities(self, velocities):
619 """Set the momenta by specifying the velocities."""
620 self.set_momenta(self.get_masses()[:, np.newaxis] * velocities)
622 def get_momenta(self):
623 """Get array of momenta."""
624 if 'momenta' in self.arrays:
625 return self.arrays['momenta'].copy()
626 else:
627 return np.zeros((len(self), 3))
629 def set_masses(self, masses='defaults'):
630 """Set atomic masses in atomic mass units.
632 The array masses should contain a list of masses. In case
633 the masses argument is not given or for those elements of the
634 masses list that are None, standard values are set."""
636 if isinstance(masses, str):
637 if masses == 'defaults':
638 masses = atomic_masses[self.arrays['numbers']]
639 elif masses == 'most_common':
640 masses = atomic_masses_common[self.arrays['numbers']]
641 elif masses is None:
642 pass
643 elif not isinstance(masses, np.ndarray):
644 masses = list(masses)
645 for i, mass in enumerate(masses):
646 if mass is None:
647 masses[i] = atomic_masses[self.numbers[i]]
648 self.set_array('masses', masses, float, ())
650 def get_masses(self):
651 """Get array of masses in atomic mass units."""
652 if 'masses' in self.arrays:
653 return self.arrays['masses'].copy()
654 else:
655 return atomic_masses[self.arrays['numbers']]
657 def set_initial_magnetic_moments(self, magmoms=None):
658 """Set the initial magnetic moments.
660 Use either one or three numbers for every atom (collinear
661 or non-collinear spins)."""
663 if magmoms is None:
664 self.set_array('initial_magmoms', None)
665 else:
666 magmoms = np.asarray(magmoms)
667 self.set_array('initial_magmoms', magmoms, float,
668 magmoms.shape[1:])
670 def get_initial_magnetic_moments(self):
671 """Get array of initial magnetic moments."""
672 if 'initial_magmoms' in self.arrays:
673 return self.arrays['initial_magmoms'].copy()
674 else:
675 return np.zeros(len(self))
677 def get_magnetic_moments(self):
678 """Get calculated local magnetic moments."""
679 if self._calc is None:
680 raise RuntimeError('Atoms object has no calculator.')
681 return self._calc.get_magnetic_moments(self)
683 def get_magnetic_moment(self):
684 """Get calculated total magnetic moment."""
685 if self._calc is None:
686 raise RuntimeError('Atoms object has no calculator.')
687 return self._calc.get_magnetic_moment(self)
689 def set_initial_charges(self, charges=None):
690 """Set the initial charges."""
692 if charges is None:
693 self.set_array('initial_charges', None)
694 else:
695 self.set_array('initial_charges', charges, float, ())
697 def get_initial_charges(self):
698 """Get array of initial charges."""
699 if 'initial_charges' in self.arrays:
700 return self.arrays['initial_charges'].copy()
701 else:
702 return np.zeros(len(self))
704 def get_charges(self):
705 """Get calculated charges."""
706 if self._calc is None:
707 raise RuntimeError('Atoms object has no calculator.')
708 try:
709 return self._calc.get_charges(self)
710 except AttributeError:
711 from ase.calculators.calculator import PropertyNotImplementedError
712 raise PropertyNotImplementedError
714 def set_positions(self, newpositions, apply_constraint=True):
715 """Set positions, honoring any constraints. To ignore constraints,
716 use *apply_constraint=False*."""
717 if self.constraints and apply_constraint:
718 newpositions = np.array(newpositions, float)
719 for constraint in self.constraints:
720 constraint.adjust_positions(self, newpositions)
722 self.set_array('positions', newpositions, shape=(3,))
724 def get_positions(self, wrap=False, **wrap_kw):
725 """Get array of positions.
727 Parameters:
729 wrap: bool
730 wrap atoms back to the cell before returning positions
731 wrap_kw: (keyword=value) pairs
732 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
733 see :func:`ase.geometry.wrap_positions`
734 """
735 from ase.geometry import wrap_positions
736 if wrap:
737 if 'pbc' not in wrap_kw:
738 wrap_kw['pbc'] = self.pbc
739 return wrap_positions(self.positions, self.cell, **wrap_kw)
740 else:
741 return self.arrays['positions'].copy()
743 def get_potential_energy(self, force_consistent=False,
744 apply_constraint=True):
745 """Calculate potential energy.
747 Ask the attached calculator to calculate the potential energy and
748 apply constraints. Use *apply_constraint=False* to get the raw
749 forces.
751 When supported by the calculator, either the energy extrapolated
752 to zero Kelvin or the energy consistent with the forces (the free
753 energy) can be returned.
754 """
755 if self._calc is None:
756 raise RuntimeError('Atoms object has no calculator.')
757 if force_consistent:
758 energy = self._calc.get_potential_energy(
759 self, force_consistent=force_consistent)
760 else:
761 energy = self._calc.get_potential_energy(self)
762 if apply_constraint:
763 for constraint in self.constraints:
764 if hasattr(constraint, 'adjust_potential_energy'):
765 energy += constraint.adjust_potential_energy(self)
766 return energy
768 def get_properties(self, properties):
769 """This method is experimental; currently for internal use."""
770 # XXX Something about constraints.
771 if self._calc is None:
772 raise RuntimeError('Atoms object has no calculator.')
773 return self._calc.calculate_properties(self, properties)
775 def get_potential_energies(self):
776 """Calculate the potential energies of all the atoms.
778 Only available with calculators supporting per-atom energies
779 (e.g. classical potentials).
780 """
781 if self._calc is None:
782 raise RuntimeError('Atoms object has no calculator.')
783 return self._calc.get_potential_energies(self)
785 def get_kinetic_energy(self):
786 """Get the kinetic energy."""
787 momenta = self.arrays.get('momenta')
788 if momenta is None:
789 return 0.0
790 return 0.5 * np.vdot(momenta, self.get_velocities())
792 def get_velocities(self):
793 """Get array of velocities."""
794 momenta = self.get_momenta()
795 masses = self.get_masses()
796 return momenta / masses[:, np.newaxis]
798 def get_total_energy(self):
799 """Get the total energy - potential plus kinetic energy."""
800 return self.get_potential_energy() + self.get_kinetic_energy()
802 def get_forces(self, apply_constraint=True, md=False):
803 """Calculate atomic forces.
805 Ask the attached calculator to calculate the forces and apply
806 constraints. Use *apply_constraint=False* to get the raw
807 forces.
809 For molecular dynamics (md=True) we don't apply the constraint
810 to the forces but to the momenta. When holonomic constraints for
811 rigid linear triatomic molecules are present, ask the constraints
812 to redistribute the forces within each triple defined in the
813 constraints (required for molecular dynamics with this type of
814 constraints)."""
816 if self._calc is None:
817 raise RuntimeError('Atoms object has no calculator.')
818 forces = self._calc.get_forces(self)
820 if apply_constraint:
821 # We need a special md flag here because for MD we want
822 # to skip real constraints but include special "constraints"
823 # Like Hookean.
824 for constraint in self.constraints:
825 if md and hasattr(constraint, 'redistribute_forces_md'):
826 constraint.redistribute_forces_md(self, forces)
827 if not md or hasattr(constraint, 'adjust_potential_energy'):
828 constraint.adjust_forces(self, forces)
829 return forces
831 # Informs calculators (e.g. Asap) that ideal gas contribution is added here.
832 _ase_handles_dynamic_stress = True
834 def get_stress(self, voigt=True, apply_constraint=True,
835 include_ideal_gas=False):
836 """Calculate stress tensor.
838 Returns an array of the six independent components of the
839 symmetric stress tensor, in the traditional Voigt order
840 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt
841 order.
843 The ideal gas contribution to the stresses is added if the
844 atoms have momenta and ``include_ideal_gas`` is set to True.
845 """
847 if self._calc is None:
848 raise RuntimeError('Atoms object has no calculator.')
850 stress = self._calc.get_stress(self)
851 shape = stress.shape
853 if shape == (3, 3):
854 # Convert to the Voigt form before possibly applying
855 # constraints and adding the dynamic part of the stress
856 # (the "ideal gas contribution").
857 stress = full_3x3_to_voigt_6_stress(stress)
858 else:
859 assert shape == (6,)
861 if apply_constraint:
862 for constraint in self.constraints:
863 if hasattr(constraint, 'adjust_stress'):
864 constraint.adjust_stress(self, stress)
866 # Add ideal gas contribution, if applicable
867 if include_ideal_gas and self.has('momenta'):
868 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
869 p = self.get_momenta()
870 masses = self.get_masses()
871 invmass = 1.0 / masses
872 invvol = 1.0 / self.get_volume()
873 for alpha in range(3):
874 for beta in range(alpha, 3):
875 stress[stresscomp[alpha, beta]] -= (
876 p[:, alpha] * p[:, beta] * invmass).sum() * invvol
878 if voigt:
879 return stress
880 else:
881 return voigt_6_to_full_3x3_stress(stress)
883 def get_stresses(self, include_ideal_gas=False, voigt=True):
884 """Calculate the stress-tensor of all the atoms.
886 Only available with calculators supporting per-atom energies and
887 stresses (e.g. classical potentials). Even for such calculators
888 there is a certain arbitrariness in defining per-atom stresses.
890 The ideal gas contribution to the stresses is added if the
891 atoms have momenta and ``include_ideal_gas`` is set to True.
892 """
893 if self._calc is None:
894 raise RuntimeError('Atoms object has no calculator.')
895 stresses = self._calc.get_stresses(self)
897 # make sure `stresses` are in voigt form
898 if np.shape(stresses)[1:] == (3, 3):
899 stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses]
900 stresses = np.array(stresses_voigt)
902 # REMARK: The ideal gas contribution is intensive, i.e., the volume
903 # is divided out. We currently don't check if `stresses` are intensive
904 # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`.
905 # It might be good to check this here, but adds computational overhead.
907 if include_ideal_gas and self.has('momenta'):
908 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
909 if hasattr(self._calc, 'get_atomic_volumes'):
910 invvol = 1.0 / self._calc.get_atomic_volumes()
911 else:
912 invvol = self.get_global_number_of_atoms() / self.get_volume()
913 p = self.get_momenta()
914 invmass = 1.0 / self.get_masses()
915 for alpha in range(3):
916 for beta in range(alpha, 3):
917 stresses[:, stresscomp[alpha, beta]] -= (
918 p[:, alpha] * p[:, beta] * invmass * invvol)
919 if voigt:
920 return stresses
921 else:
922 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
923 return np.array(stresses_3x3)
925 def get_dipole_moment(self):
926 """Calculate the electric dipole moment for the atoms object.
928 Only available for calculators which has a get_dipole_moment()
929 method."""
931 if self._calc is None:
932 raise RuntimeError('Atoms object has no calculator.')
933 return self._calc.get_dipole_moment(self)
935 def copy(self):
936 """Return a copy."""
937 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
938 celldisp=self._celldisp.copy())
940 atoms.arrays = {}
941 for name, a in self.arrays.items():
942 atoms.arrays[name] = a.copy()
943 atoms.constraints = copy.deepcopy(self.constraints)
944 return atoms
946 def todict(self):
947 """For basic JSON (non-database) support."""
948 d = dict(self.arrays)
949 d['cell'] = np.asarray(self.cell)
950 d['pbc'] = self.pbc
951 if self._celldisp.any():
952 d['celldisp'] = self._celldisp
953 if self.constraints:
954 d['constraints'] = self.constraints
955 if self.info:
956 d['info'] = self.info
957 # Calculator... trouble.
958 return d
960 @classmethod
961 def fromdict(cls, dct):
962 """Rebuild atoms object from dictionary representation (todict)."""
963 dct = dct.copy()
964 kw = {name: dct.pop(name)
965 for name in ['numbers', 'positions', 'cell', 'pbc']}
966 constraints = dct.pop('constraints', None)
967 if constraints:
968 from ase.constraints import dict2constraint
969 constraints = [dict2constraint(d) for d in constraints]
971 info = dct.pop('info', None)
973 atoms = cls(constraint=constraints,
974 celldisp=dct.pop('celldisp', None),
975 info=info, **kw)
976 natoms = len(atoms)
978 # Some arrays are named differently from the atoms __init__ keywords.
979 # Also, there may be custom arrays. Hence we set them directly:
980 for name, arr in dct.items():
981 assert len(arr) == natoms, name
982 assert isinstance(arr, np.ndarray)
983 atoms.arrays[name] = arr
984 return atoms
986 def __len__(self):
987 return len(self.arrays['positions'])
989 @deprecated(
990 "Please use len(self) or, if your atoms are distributed, "
991 "self.get_global_number_of_atoms.",
992 category=FutureWarning,
993 )
994 def get_number_of_atoms(self):
995 """
996 .. deprecated:: 3.18.1
997 You probably want ``len(atoms)``. Or if your atoms are distributed,
998 use (and see) :func:`get_global_number_of_atoms()`.
999 """
1000 return len(self)
1002 def get_global_number_of_atoms(self):
1003 """Returns the global number of atoms in a distributed-atoms parallel
1004 simulation.
1006 DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING!
1008 Equivalent to len(atoms) in the standard ASE Atoms class. You should
1009 normally use len(atoms) instead. This function's only purpose is to
1010 make compatibility between ASE and Asap easier to maintain by having a
1011 few places in ASE use this function instead. It is typically only
1012 when counting the global number of degrees of freedom or in similar
1013 situations.
1014 """
1015 return len(self)
1017 def __repr__(self):
1018 tokens = []
1020 N = len(self)
1021 if N <= 60:
1022 symbols = self.get_chemical_formula('reduce')
1023 else:
1024 symbols = self.get_chemical_formula('hill')
1025 tokens.append(f"symbols='{symbols}'")
1027 if self.pbc.any() and not self.pbc.all():
1028 tokens.append(f'pbc={self.pbc.tolist()}')
1029 else:
1030 tokens.append(f'pbc={self.pbc[0]}')
1032 cell = self.cell
1033 if cell:
1034 if cell.orthorhombic:
1035 cell = cell.lengths().tolist()
1036 else:
1037 cell = cell.tolist()
1038 tokens.append(f'cell={cell}')
1040 for name in sorted(self.arrays):
1041 if name in ['numbers', 'positions']:
1042 continue
1043 tokens.append(f'{name}=...')
1045 if self.constraints:
1046 if len(self.constraints) == 1:
1047 constraint = self.constraints[0]
1048 else:
1049 constraint = self.constraints
1050 tokens.append(f'constraint={constraint!r}')
1052 if self._calc is not None:
1053 tokens.append('calculator={}(...)'
1054 .format(self._calc.__class__.__name__))
1056 return '{}({})'.format(self.__class__.__name__, ', '.join(tokens))
1058 def __add__(self, other):
1059 atoms = self.copy()
1060 atoms += other
1061 return atoms
1063 def extend(self, other):
1064 """Extend atoms object by appending atoms from *other*."""
1065 if isinstance(other, Atom):
1066 other = self.__class__([other])
1068 n1 = len(self)
1069 n2 = len(other)
1071 for name, a1 in self.arrays.items():
1072 a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype)
1073 a[:n1] = a1
1074 if name == 'masses':
1075 a2 = other.get_masses()
1076 else:
1077 a2 = other.arrays.get(name)
1078 if a2 is not None:
1079 a[n1:] = a2
1080 self.arrays[name] = a
1082 for name, a2 in other.arrays.items():
1083 if name in self.arrays:
1084 continue
1085 a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype)
1086 a[n1:] = a2
1087 if name == 'masses':
1088 a[:n1] = self.get_masses()[:n1]
1089 else:
1090 a[:n1] = 0
1092 self.set_array(name, a)
1094 def __iadd__(self, other):
1095 self.extend(other)
1096 return self
1098 def append(self, atom):
1099 """Append atom to end."""
1100 self.extend(self.__class__([atom]))
1102 def __iter__(self):
1103 for i in range(len(self)):
1104 yield self[i]
1106 def __getitem__(self, i):
1107 """Return a subset of the atoms.
1109 i -- scalar integer, list of integers, or slice object
1110 describing which atoms to return.
1112 If i is a scalar, return an Atom object. If i is a list or a
1113 slice, return an Atoms object with the same cell, pbc, and
1114 other associated info as the original Atoms object. The
1115 indices of the constraints will be shuffled so that they match
1116 the indexing in the subset returned.
1118 """
1120 if isinstance(i, numbers.Integral):
1121 natoms = len(self)
1122 if i < -natoms or i >= natoms:
1123 raise IndexError('Index out of range.')
1125 return Atom(atoms=self, index=i)
1126 elif not isinstance(i, slice):
1127 i = np.array(i)
1128 if len(i) == 0:
1129 i = np.array([], dtype=int)
1130 # if i is a mask
1131 if i.dtype == bool:
1132 if len(i) != len(self):
1133 raise IndexError('Length of mask {} must equal '
1134 'number of atoms {}'
1135 .format(len(i), len(self)))
1136 i = np.arange(len(self))[i]
1138 import copy
1140 conadd = []
1141 # Constraints need to be deepcopied, but only the relevant ones.
1142 for con in copy.deepcopy(self.constraints):
1143 try:
1144 con.index_shuffle(self, i)
1145 except (IndexError, NotImplementedError):
1146 pass
1147 else:
1148 conadd.append(con)
1150 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1151 # should be communicated to the slice as well
1152 celldisp=self._celldisp)
1153 # TODO: Do we need to shuffle indices in adsorbate_info too?
1155 atoms.arrays = {}
1156 for name, a in self.arrays.items():
1157 atoms.arrays[name] = a[i].copy()
1159 atoms.constraints = conadd
1160 return atoms
1162 def __delitem__(self, i):
1163 from ase.constraints import FixAtoms
1164 for c in self._constraints:
1165 if not isinstance(c, FixAtoms):
1166 raise RuntimeError('Remove constraint using set_constraint() '
1167 'before deleting atoms.')
1169 if isinstance(i, list) and len(i) > 0:
1170 # Make sure a list of booleans will work correctly and not be
1171 # interpreted at 0 and 1 indices.
1172 i = np.array(i)
1174 if len(self._constraints) > 0:
1175 n = len(self)
1176 i = np.arange(n)[i]
1177 if isinstance(i, int):
1178 i = [i]
1179 constraints = []
1180 for c in self._constraints:
1181 c = c.delete_atoms(i, n)
1182 if c is not None:
1183 constraints.append(c)
1184 self.constraints = constraints
1186 mask = np.ones(len(self), bool)
1187 mask[i] = False
1188 for name, a in self.arrays.items():
1189 self.arrays[name] = a[mask]
1191 def pop(self, i=-1):
1192 """Remove and return atom at index *i* (default last)."""
1193 atom = self[i]
1194 atom.cut_reference_to_atoms()
1195 del self[i]
1196 return atom
1198 def __imul__(self, m):
1199 """In-place repeat of atoms."""
1200 if isinstance(m, int):
1201 m = (m, m, m)
1203 for x, vec in zip(m, self.cell):
1204 if x != 1 and not vec.any():
1205 raise ValueError('Cannot repeat along undefined lattice '
1206 'vector')
1208 M = np.prod(m)
1209 n = len(self)
1211 for name, a in self.arrays.items():
1212 self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1))
1214 positions = self.arrays['positions']
1215 i0 = 0
1216 for m0 in range(m[0]):
1217 for m1 in range(m[1]):
1218 for m2 in range(m[2]):
1219 i1 = i0 + n
1220 positions[i0:i1] += np.dot((m0, m1, m2), self.cell)
1221 i0 = i1
1223 if self.constraints is not None:
1224 self.constraints = [c.repeat(m, n) for c in self.constraints]
1226 self.cell = np.array([m[c] * self.cell[c] for c in range(3)])
1228 return self
1230 def repeat(self, rep):
1231 """Create new repeated atoms object.
1233 The *rep* argument should be a sequence of three positive
1234 integers like *(2,3,1)* or a single integer (*r*) equivalent
1235 to *(r,r,r)*."""
1237 atoms = self.copy()
1238 atoms *= rep
1239 return atoms
1241 def __mul__(self, rep):
1242 return self.repeat(rep)
1244 def translate(self, displacement):
1245 """Translate atomic positions.
1247 The displacement argument can be a float an xyz vector or an
1248 nx3 array (where n is the number of atoms)."""
1250 self.arrays['positions'] += np.array(displacement)
1252 def center(self, vacuum=None, axis=(0, 1, 2), about=None):
1253 """Center atoms in unit cell.
1255 Centers the atoms in the unit cell, so there is the same
1256 amount of vacuum on all sides.
1258 vacuum: float (default: None)
1259 If specified adjust the amount of vacuum when centering.
1260 If vacuum=10.0 there will thus be 10 Angstrom of vacuum
1261 on each side.
1262 axis: int or sequence of ints
1263 Axis or axes to act on. Default: Act on all axes.
1264 about: float or array (default: None)
1265 If specified, center the atoms about <about>.
1266 I.e., about=(0., 0., 0.) (or just "about=0.", interpreted
1267 identically), to center about the origin.
1268 """
1270 # Find the orientations of the faces of the unit cell
1271 cell = self.cell.complete()
1272 dirs = np.zeros_like(cell)
1274 lengths = cell.lengths()
1275 for i in range(3):
1276 dirs[i] = np.cross(cell[i - 1], cell[i - 2])
1277 dirs[i] /= np.linalg.norm(dirs[i])
1278 if dirs[i] @ cell[i] < 0.0:
1279 dirs[i] *= -1
1281 if isinstance(axis, int):
1282 axes = (axis,)
1283 else:
1284 axes = axis
1286 # Now, decide how much each basis vector should be made longer
1287 pos = self.positions
1288 longer = np.zeros(3)
1289 shift = np.zeros(3)
1290 for i in axes:
1291 if len(pos):
1292 scalarprod = pos @ dirs[i]
1293 p0 = scalarprod.min()
1294 p1 = scalarprod.max()
1295 else:
1296 p0 = 0
1297 p1 = 0
1298 height = cell[i] @ dirs[i]
1299 if vacuum is not None:
1300 lng = (p1 - p0 + 2 * vacuum) - height
1301 else:
1302 lng = 0.0 # Do not change unit cell size!
1303 top = lng + height - p1
1304 shf = 0.5 * (top - p0)
1305 cosphi = cell[i] @ dirs[i] / lengths[i]
1306 longer[i] = lng / cosphi
1307 shift[i] = shf / cosphi
1309 # Now, do it!
1310 translation = np.zeros(3)
1311 for i in axes:
1312 nowlen = lengths[i]
1313 if vacuum is not None:
1314 self.cell[i] = cell[i] * (1 + longer[i] / nowlen)
1315 translation += shift[i] * cell[i] / nowlen
1317 # We calculated translations using the completed cell,
1318 # so directions without cell vectors will have been centered
1319 # along a "fake" vector of length 1.
1320 # Therefore, we adjust by -0.5:
1321 if not any(self.cell[i]):
1322 translation[i] -= 0.5
1324 # Optionally, translate to center about a point in space.
1325 if about is not None:
1326 for n, vector in enumerate(self.cell):
1327 if n in axes:
1328 translation -= vector / 2.0
1329 translation[n] += about[n]
1331 self.positions += translation
1333 def get_center_of_mass(self, scaled=False, indices=None):
1334 """Get the center of mass.
1336 Parameters
1337 ----------
1338 scaled : bool
1339 If True, the center of mass in scaled coordinates is returned.
1340 indices : list | slice | str, default: None
1341 If specified, the center of mass of a subset of atoms is returned.
1342 """
1343 if indices is None:
1344 indices = slice(None)
1345 elif isinstance(indices, str):
1346 indices = string2index(indices)
1348 masses = self.get_masses()[indices]
1349 com = masses @ self.positions[indices] / masses.sum()
1350 if scaled:
1351 return self.cell.scaled_positions(com)
1352 return com # Cartesian coordinates
1354 def set_center_of_mass(self, com, scaled=False):
1355 """Set the center of mass.
1357 If scaled=True the center of mass is expected in scaled coordinates.
1358 Constraints are considered for scaled=False.
1359 """
1360 old_com = self.get_center_of_mass(scaled=scaled)
1361 difference = com - old_com
1362 if scaled:
1363 self.set_scaled_positions(self.get_scaled_positions() + difference)
1364 else:
1365 self.set_positions(self.get_positions() + difference)
1367 def get_moments_of_inertia(self, vectors=False):
1368 """Get the moments of inertia along the principal axes.
1370 The three principal moments of inertia are computed from the
1371 eigenvalues of the symmetric inertial tensor. Periodic boundary
1372 conditions are ignored. Units of the moments of inertia are
1373 amu*angstrom**2.
1374 """
1375 com = self.get_center_of_mass()
1376 positions = self.get_positions()
1377 positions -= com # translate center of mass to origin
1378 masses = self.get_masses()
1380 # Initialize elements of the inertial tensor
1381 I11 = I22 = I33 = I12 = I13 = I23 = 0.0
1382 for i in range(len(self)):
1383 x, y, z = positions[i]
1384 m = masses[i]
1386 I11 += m * (y ** 2 + z ** 2)
1387 I22 += m * (x ** 2 + z ** 2)
1388 I33 += m * (x ** 2 + y ** 2)
1389 I12 += -m * x * y
1390 I13 += -m * x * z
1391 I23 += -m * y * z
1393 Itensor = np.array([[I11, I12, I13],
1394 [I12, I22, I23],
1395 [I13, I23, I33]])
1397 evals, evecs = np.linalg.eigh(Itensor)
1398 if vectors:
1399 return evals, evecs.transpose()
1400 else:
1401 return evals
1403 def get_angular_momentum(self):
1404 """Get total angular momentum with respect to the center of mass."""
1405 com = self.get_center_of_mass()
1406 positions = self.get_positions()
1407 positions -= com # translate center of mass to origin
1408 return np.cross(positions, self.get_momenta()).sum(0)
1410 def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False):
1411 """Rotate atoms based on a vector and an angle, or two vectors.
1413 Parameters:
1415 a = None:
1416 Angle that the atoms is rotated around the vector 'v'. 'a'
1417 can also be a vector and then 'a' is rotated
1418 into 'v'.
1420 v:
1421 Vector to rotate the atoms around. Vectors can be given as
1422 strings: 'x', '-x', 'y', ... .
1424 center = (0, 0, 0):
1425 The center is kept fixed under the rotation. Use 'COM' to fix
1426 the center of mass, 'COP' to fix the center of positions or
1427 'COU' to fix the center of cell.
1429 rotate_cell = False:
1430 If true the cell is also rotated.
1432 Examples:
1434 Rotate 90 degrees around the z-axis, so that the x-axis is
1435 rotated into the y-axis:
1437 >>> atoms = Atoms()
1438 >>> atoms.rotate(90, 'z')
1439 >>> atoms.rotate(90, (0, 0, 1))
1440 >>> atoms.rotate(-90, '-z')
1441 >>> atoms.rotate('x', 'y')
1442 >>> atoms.rotate((1, 0, 0), (0, 1, 0))
1443 """
1445 if not isinstance(a, numbers.Real):
1446 a, v = v, a
1448 norm = np.linalg.norm
1449 v = string2vector(v)
1451 normv = norm(v)
1453 if normv == 0.0:
1454 raise ZeroDivisionError('Cannot rotate: norm(v) == 0')
1456 if isinstance(a, numbers.Real):
1457 a *= pi / 180
1458 v /= normv
1459 c = cos(a)
1460 s = sin(a)
1461 else:
1462 v2 = string2vector(a)
1463 v /= normv
1464 normv2 = np.linalg.norm(v2)
1465 if normv2 == 0:
1466 raise ZeroDivisionError('Cannot rotate: norm(a) == 0')
1467 v2 /= norm(v2)
1468 c = np.dot(v, v2)
1469 v = np.cross(v, v2)
1470 s = norm(v)
1471 # In case *v* and *a* are parallel, np.cross(v, v2) vanish
1472 # and can't be used as a rotation axis. However, in this
1473 # case any rotation axis perpendicular to v2 will do.
1474 eps = 1e-7
1475 if s < eps:
1476 v = np.cross((0, 0, 1), v2)
1477 if norm(v) < eps:
1478 v = np.cross((1, 0, 0), v2)
1479 assert norm(v) >= eps
1480 elif s > 0:
1481 v /= s
1483 center = self._centering_as_array(center)
1485 p = self.arrays['positions'] - center
1486 self.arrays['positions'][:] = (c * p -
1487 np.cross(p, s * v) +
1488 np.outer(np.dot(p, v), (1.0 - c) * v) +
1489 center)
1490 if rotate_cell:
1491 rotcell = self.get_cell()
1492 rotcell[:] = (c * rotcell -
1493 np.cross(rotcell, s * v) +
1494 np.outer(np.dot(rotcell, v), (1.0 - c) * v))
1495 self.set_cell(rotcell)
1497 def _centering_as_array(self, center):
1498 if isinstance(center, str):
1499 if center.lower() == 'com':
1500 center = self.get_center_of_mass()
1501 elif center.lower() == 'cop':
1502 center = self.get_positions().mean(axis=0)
1503 elif center.lower() == 'cou':
1504 center = self.get_cell().sum(axis=0) / 2
1505 else:
1506 raise ValueError('Cannot interpret center')
1507 else:
1508 center = np.array(center, float)
1509 return center
1511 def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)):
1512 """Rotate atoms via Euler angles (in degrees).
1514 See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation.
1516 Parameters:
1518 center :
1519 The point to rotate about. A sequence of length 3 with the
1520 coordinates, or 'COM' to select the center of mass, 'COP' to
1521 select center of positions or 'COU' to select center of cell.
1522 phi :
1523 The 1st rotation angle around the z axis.
1524 theta :
1525 Rotation around the x axis.
1526 psi :
1527 2nd rotation around the z axis.
1529 """
1530 center = self._centering_as_array(center)
1532 phi *= pi / 180
1533 theta *= pi / 180
1534 psi *= pi / 180
1536 # First move the molecule to the origin In contrast to MATLAB,
1537 # numpy broadcasts the smaller array to the larger row-wise,
1538 # so there is no need to play with the Kronecker product.
1539 rcoords = self.positions - center
1540 # First Euler rotation about z in matrix form
1541 D = np.array(((cos(phi), sin(phi), 0.),
1542 (-sin(phi), cos(phi), 0.),
1543 (0., 0., 1.)))
1544 # Second Euler rotation about x:
1545 C = np.array(((1., 0., 0.),
1546 (0., cos(theta), sin(theta)),
1547 (0., -sin(theta), cos(theta))))
1548 # Third Euler rotation, 2nd rotation about z:
1549 B = np.array(((cos(psi), sin(psi), 0.),
1550 (-sin(psi), cos(psi), 0.),
1551 (0., 0., 1.)))
1552 # Total Euler rotation
1553 A = np.dot(B, np.dot(C, D))
1554 # Do the rotation
1555 rcoords = np.dot(A, np.transpose(rcoords))
1556 # Move back to the rotation point
1557 self.positions = np.transpose(rcoords) + center
1559 def get_dihedral(self, a0, a1, a2, a3, mic=False):
1560 """Calculate dihedral angle.
1562 Calculate dihedral angle (in degrees) between the vectors a0->a1
1563 and a2->a3.
1565 Use mic=True to use the Minimum Image Convention and calculate the
1566 angle across periodic boundaries.
1567 """
1568 return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0]
1570 def get_dihedrals(self, indices, mic=False):
1571 """Calculate dihedral angles.
1573 Calculate dihedral angles (in degrees) between the list of vectors
1574 a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices.
1576 Use mic=True to use the Minimum Image Convention and calculate the
1577 angles across periodic boundaries.
1578 """
1579 from ase.geometry import get_dihedrals
1581 indices = np.array(indices)
1582 assert indices.shape[1] == 4
1584 a0s = self.positions[indices[:, 0]]
1585 a1s = self.positions[indices[:, 1]]
1586 a2s = self.positions[indices[:, 2]]
1587 a3s = self.positions[indices[:, 3]]
1589 # vectors 0->1, 1->2, 2->3
1590 v0 = a1s - a0s
1591 v1 = a2s - a1s
1592 v2 = a3s - a2s
1594 cell = None
1595 pbc = None
1597 if mic:
1598 cell = self.cell
1599 pbc = self.pbc
1601 return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc)
1603 def _masked_rotate(self, center, axis, diff, mask):
1604 # do rotation of subgroup by copying it to temporary atoms object
1605 # and then rotating that
1606 #
1607 # recursive object definition might not be the most elegant thing,
1608 # more generally useful might be a rotation function with a mask?
1609 group = self.__class__()
1610 for i in range(len(self)):
1611 if mask[i]:
1612 group += self[i]
1613 group.translate(-center)
1614 group.rotate(diff * 180 / pi, axis)
1615 group.translate(center)
1616 # set positions in original atoms object
1617 j = 0
1618 for i in range(len(self)):
1619 if mask[i]:
1620 self.positions[i] = group[j].position
1621 j += 1
1623 def set_dihedral(self, a1, a2, a3, a4, angle,
1624 mask=None, indices=None):
1625 """Set the dihedral angle (degrees) between vectors a1->a2 and
1626 a3->a4 by changing the atom indexed by a4.
1628 If mask is not None, all the atoms described in mask
1629 (read: the entire subgroup) are moved. Alternatively to the mask,
1630 the indices of the atoms to be rotated can be supplied. If both
1631 *mask* and *indices* are given, *indices* overwrites *mask*.
1633 **Important**: If *mask* or *indices* is given and does not contain
1634 *a4*, *a4* will NOT be moved. In most cases you therefore want
1635 to include *a4* in *mask*/*indices*.
1637 Example: the following defines a very crude
1638 ethane-like molecule and twists one half of it by 30 degrees.
1640 >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0],
1641 ... [1, 0, 0], [2, 1, 0], [2, -1, 0]])
1642 >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1])
1643 """
1645 angle *= pi / 180
1647 # if not provided, set mask to the last atom in the
1648 # dihedral description
1649 if mask is None and indices is None:
1650 mask = np.zeros(len(self))
1651 mask[a4] = 1
1652 elif indices is not None:
1653 mask = [index in indices for index in range(len(self))]
1655 # compute necessary in dihedral change, from current value
1656 current = self.get_dihedral(a1, a2, a3, a4) * pi / 180
1657 diff = angle - current
1658 axis = self.positions[a3] - self.positions[a2]
1659 center = self.positions[a3]
1660 self._masked_rotate(center, axis, diff, mask)
1662 def rotate_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None):
1663 """Rotate dihedral angle.
1665 Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a
1666 predefined dihedral angle, starting from its current configuration.
1667 """
1668 start = self.get_dihedral(a1, a2, a3, a4)
1669 self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices)
1671 def get_angle(self, a1, a2, a3, mic=False):
1672 """Get angle formed by three atoms.
1674 Calculate angle in degrees between the vectors a2->a1 and
1675 a2->a3.
1677 Use mic=True to use the Minimum Image Convention and calculate the
1678 angle across periodic boundaries.
1679 """
1680 return self.get_angles([[a1, a2, a3]], mic=mic)[0]
1682 def get_angles(self, indices, mic=False):
1683 """Get angle formed by three atoms for multiple groupings.
1685 Calculate angle in degrees between vectors between atoms a2->a1
1686 and a2->a3, where a1, a2, and a3 are in each row of indices.
1688 Use mic=True to use the Minimum Image Convention and calculate
1689 the angle across periodic boundaries.
1690 """
1691 from ase.geometry import get_angles
1693 indices = np.array(indices)
1694 assert indices.shape[1] == 3
1696 a1s = self.positions[indices[:, 0]]
1697 a2s = self.positions[indices[:, 1]]
1698 a3s = self.positions[indices[:, 2]]
1700 v12 = a1s - a2s
1701 v32 = a3s - a2s
1703 cell = None
1704 pbc = None
1706 if mic:
1707 cell = self.cell
1708 pbc = self.pbc
1710 return get_angles(v12, v32, cell=cell, pbc=pbc)
1712 def set_angle(self, a1, a2=None, a3=None, angle=None, mask=None,
1713 indices=None, add=False):
1714 """Set angle (in degrees) formed by three atoms.
1716 Sets the angle between vectors *a2*->*a1* and *a2*->*a3*.
1718 If *add* is `True`, the angle will be changed by the value given.
1720 Same usage as in :meth:`ase.Atoms.set_dihedral`.
1721 If *mask* and *indices*
1722 are given, *indices* overwrites *mask*. If *mask* and *indices*
1723 are not set, only *a3* is moved."""
1725 if any(a is None for a in [a2, a3, angle]):
1726 raise ValueError('a2, a3, and angle must not be None')
1728 # If not provided, set mask to the last atom in the angle description
1729 if mask is None and indices is None:
1730 mask = np.zeros(len(self))
1731 mask[a3] = 1
1732 elif indices is not None:
1733 mask = [index in indices for index in range(len(self))]
1735 if add:
1736 diff = angle
1737 else:
1738 # Compute necessary in angle change, from current value
1739 diff = angle - self.get_angle(a1, a2, a3)
1741 diff *= pi / 180
1742 # Do rotation of subgroup by copying it to temporary atoms object and
1743 # then rotating that
1744 v10 = self.positions[a1] - self.positions[a2]
1745 v12 = self.positions[a3] - self.positions[a2]
1746 v10 /= np.linalg.norm(v10)
1747 v12 /= np.linalg.norm(v12)
1748 axis = np.cross(v10, v12)
1749 center = self.positions[a2]
1750 self._masked_rotate(center, axis, diff, mask)
1752 def rattle(self, stdev=0.001, seed=None, rng=None):
1753 """Randomly displace atoms.
1755 This method adds random displacements to the atomic positions,
1756 taking a possible constraint into account. The random numbers are
1757 drawn from a normal distribution of standard deviation stdev.
1759 By default, the random number generator always uses the same seed (42)
1760 for repeatability. You can provide your own seed (an integer), or if you
1761 want the randomness to be different each time you run a script, then
1762 provide `rng=numpy.random`. For a parallel calculation, it is important
1763 to use the same seed on all processors! """
1765 if seed is not None and rng is not None:
1766 raise ValueError('Please do not provide both seed and rng.')
1768 if rng is None:
1769 if seed is None:
1770 seed = 42
1771 rng = np.random.RandomState(seed)
1772 positions = self.arrays['positions']
1773 self.set_positions(positions +
1774 rng.normal(scale=stdev, size=positions.shape))
1776 def get_distance(self, a0, a1, mic=False, vector=False):
1777 """Return distance between two atoms.
1779 Use mic=True to use the Minimum Image Convention.
1780 vector=True gives the distance vector (from a0 to a1).
1781 """
1782 return self.get_distances(a0, [a1], mic=mic, vector=vector)[0]
1784 def get_distances(self, a, indices, mic=False, vector=False):
1785 """Return distances of atom No.i with a list of atoms.
1787 Use mic=True to use the Minimum Image Convention.
1788 vector=True gives the distance vector (from a to self[indices]).
1789 """
1790 from ase.geometry import get_distances
1792 R = self.arrays['positions']
1793 p1 = [R[a]]
1794 p2 = R[indices]
1796 cell = None
1797 pbc = None
1799 if mic:
1800 cell = self.cell
1801 pbc = self.pbc
1803 D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc)
1805 if vector:
1806 D.shape = (-1, 3)
1807 return D
1808 else:
1809 D_len.shape = (-1,)
1810 return D_len
1812 def get_all_distances(self, mic=False, vector=False):
1813 """Return distances of all of the atoms with all of the atoms.
1815 Use mic=True to use the Minimum Image Convention.
1816 """
1817 from ase.geometry import get_distances
1819 R = self.arrays['positions']
1821 cell = None
1822 pbc = None
1824 if mic:
1825 cell = self.cell
1826 pbc = self.pbc
1828 D, D_len = get_distances(R, cell=cell, pbc=pbc)
1830 if vector:
1831 return D
1832 else:
1833 return D_len
1835 def set_distance(self, a0, a1, distance, fix=0.5, mic=False,
1836 mask=None, indices=None, add=False, factor=False):
1837 """Set the distance between two atoms.
1839 Set the distance between atoms *a0* and *a1* to *distance*.
1840 By default, the center of the two atoms will be fixed. Use
1841 *fix=0* to fix the first atom, *fix=1* to fix the second
1842 atom and *fix=0.5* (default) to fix the center of the bond.
1844 If *mask* or *indices* are set (*mask* overwrites *indices*),
1845 only the atoms defined there are moved
1846 (see :meth:`ase.Atoms.set_dihedral`).
1848 When *add* is true, the distance is changed by the value given.
1849 In combination
1850 with *factor* True, the value given is a factor scaling the distance.
1852 It is assumed that the atoms in *mask*/*indices* move together
1853 with *a1*. If *fix=1*, only *a0* will therefore be moved."""
1854 from ase.geometry import find_mic
1856 if a0 % len(self) == a1 % len(self):
1857 raise ValueError('a0 and a1 must not be the same')
1859 if add:
1860 oldDist = self.get_distance(a0, a1, mic=mic)
1861 if factor:
1862 newDist = oldDist * distance
1863 else:
1864 newDist = oldDist + distance
1865 self.set_distance(a0, a1, newDist, fix=fix, mic=mic,
1866 mask=mask, indices=indices, add=False,
1867 factor=False)
1868 return
1870 R = self.arrays['positions']
1871 D = np.array([R[a1] - R[a0]])
1873 if mic:
1874 D, D_len = find_mic(D, self.cell, self.pbc)
1875 else:
1876 D_len = np.array([np.sqrt((D**2).sum())])
1877 x = 1.0 - distance / D_len[0]
1879 if mask is None and indices is None:
1880 indices = [a0, a1]
1881 elif mask:
1882 indices = [i for i in range(len(self)) if mask[i]]
1884 for i in indices:
1885 if i == a0:
1886 R[a0] += (x * fix) * D[0]
1887 else:
1888 R[i] -= (x * (1.0 - fix)) * D[0]
1890 def get_scaled_positions(self, wrap=True):
1891 """Get positions relative to unit cell.
1893 If wrap is True, atoms outside the unit cell will be wrapped into
1894 the cell in those directions with periodic boundary conditions
1895 so that the scaled coordinates are between zero and one.
1897 If any cell vectors are zero, the corresponding coordinates
1898 are evaluated as if the cell were completed using
1899 ``cell.complete()``. This means coordinates will be Cartesian
1900 as long as the non-zero cell vectors span a Cartesian axis or
1901 plane."""
1903 fractional = self.cell.scaled_positions(self.positions)
1905 if wrap:
1906 for i, periodic in enumerate(self.pbc):
1907 if periodic:
1908 # Yes, we need to do it twice.
1909 # See the scaled_positions.py test.
1910 fractional[:, i] %= 1.0
1911 fractional[:, i] %= 1.0
1913 return fractional
1915 def set_scaled_positions(self, scaled):
1916 """Set positions relative to unit cell."""
1917 self.positions[:] = self.cell.cartesian_positions(scaled)
1919 def wrap(self, **wrap_kw):
1920 """Wrap positions to unit cell.
1922 Parameters:
1924 wrap_kw: (keyword=value) pairs
1925 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
1926 see :func:`ase.geometry.wrap_positions`
1927 """
1929 if 'pbc' not in wrap_kw:
1930 wrap_kw['pbc'] = self.pbc
1932 self.positions[:] = self.get_positions(wrap=True, **wrap_kw)
1934 def get_temperature(self):
1935 """Get the temperature in Kelvin."""
1936 ekin = self.get_kinetic_energy()
1937 return 2 * ekin / (self.get_number_of_degrees_of_freedom() * units.kB)
1939 def __eq__(self, other):
1940 """Check for identity of two atoms objects.
1942 Identity means: same positions, atomic numbers, unit cell and
1943 periodic boundary conditions."""
1944 if not isinstance(other, Atoms):
1945 return False
1946 a = self.arrays
1947 b = other.arrays
1948 return (len(self) == len(other) and
1949 (a['positions'] == b['positions']).all() and
1950 (a['numbers'] == b['numbers']).all() and
1951 (self.cell == other.cell).all() and
1952 (self.pbc == other.pbc).all())
1954 def __ne__(self, other):
1955 """Check if two atoms objects are not equal.
1957 Any differences in positions, atomic numbers, unit cell or
1958 periodic boundary condtions make atoms objects not equal.
1959 """
1960 eq = self.__eq__(other)
1961 if eq is NotImplemented:
1962 return eq
1963 else:
1964 return not eq
1966 # @deprecated('Please use atoms.cell.volume')
1967 # We kind of want to deprecate this, but the ValueError behaviour
1968 # might be desirable. Should we do this?
1969 def get_volume(self):
1970 """Get volume of unit cell."""
1971 if self.cell.rank != 3:
1972 raise ValueError(
1973 'You have {} lattice vectors: volume not defined'
1974 .format(self.cell.rank))
1975 return self.cell.volume
1977 def _get_positions(self):
1978 """Return reference to positions-array for in-place manipulations."""
1979 return self.arrays['positions']
1981 def _set_positions(self, pos):
1982 """Set positions directly, bypassing constraints."""
1983 self.arrays['positions'][:] = pos
1985 positions = property(_get_positions, _set_positions,
1986 doc='Attribute for direct ' +
1987 'manipulation of the positions.')
1989 def _get_atomic_numbers(self):
1990 """Return reference to atomic numbers for in-place
1991 manipulations."""
1992 return self.arrays['numbers']
1994 numbers = property(_get_atomic_numbers, set_atomic_numbers,
1995 doc='Attribute for direct ' +
1996 'manipulation of the atomic numbers.')
1998 @property
1999 def cell(self):
2000 """The :class:`ase.cell.Cell` for direct manipulation."""
2001 return self._cellobj
2003 @cell.setter
2004 def cell(self, cell):
2005 cell = Cell.ascell(cell)
2006 self._cellobj[:] = cell
2008 def write(self, filename, format=None, **kwargs):
2009 """Write atoms object to a file.
2011 see ase.io.write for formats.
2012 kwargs are passed to ase.io.write.
2013 """
2014 from ase.io import write
2015 write(filename, self, format, **kwargs)
2017 def iterimages(self):
2018 yield self
2020 def __ase_optimizable__(self):
2021 from ase.optimize.optimize import OptimizableAtoms
2022 return OptimizableAtoms(self)
2024 def edit(self):
2025 """Modify atoms interactively through ASE's GUI viewer.
2027 Conflicts leading to undesirable behaviour might arise
2028 when matplotlib has been pre-imported with certain
2029 incompatible backends and while trying to use the
2030 plot feature inside the interactive GUI. To circumvent,
2031 please set matplotlib.use('gtk') before calling this
2032 method.
2033 """
2034 from ase.gui.gui import GUI
2035 from ase.gui.images import Images
2036 images = Images([self])
2037 gui = GUI(images)
2038 gui.run()
2041def string2vector(v):
2042 if isinstance(v, str):
2043 if v[0] == '-':
2044 return -string2vector(v[1:])
2045 w = np.zeros(3)
2046 w['xyz'.index(v)] = 1.0
2047 return w
2048 return np.array(v, float)
2051def default(data, dflt):
2052 """Helper function for setting default values."""
2053 if data is None:
2054 return None
2055 elif isinstance(data, (list, tuple)):
2056 newdata = []
2057 allnone = True
2058 for x in data:
2059 if x is None:
2060 newdata.append(dflt)
2061 else:
2062 newdata.append(x)
2063 allnone = False
2064 if allnone:
2065 return None
2066 return newdata
2067 else:
2068 return data