Coverage for /builds/ase/ase/ase/atoms.py : 93.35%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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, sin, pi
13import numpy as np
15import ase.units as units
16from ase.atom import Atom
17from ase.cell import Cell
18from ase.stress import voigt_6_to_full_3x3_stress, full_3x3_to_voigt_6_stress
19from ase.data import atomic_masses, atomic_masses_common
20from ase.symbols import Symbols, symbols2numbers
21from ase.utils import deprecated
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 >>> d = 1.104 # N2 bondlength
107 >>> a = Atoms('N2', [(0, 0, 0), (0, 0, d)])
108 >>> a = Atoms(numbers=[7, 7], positions=[(0, 0, 0), (0, 0, d)])
109 >>> a = Atoms([Atom('N', (0, 0, 0)), Atom('N', (0, 0, d))])
111 FCC gold:
113 >>> a = 4.05 # Gold lattice constant
114 >>> b = a / 2
115 >>> fcc = Atoms('Au',
116 ... cell=[(0, b, b), (b, 0, b), (b, b, 0)],
117 ... pbc=True)
119 Hydrogen wire:
121 >>> d = 0.9 # H-H distance
122 >>> h = Atoms('H', positions=[(0, 0, 0)],
123 ... cell=(d, 0, 0),
124 ... pbc=(1, 0, 0))
125 """
127 ase_objtype = 'atoms' # For JSONability
129 def __init__(self, symbols=None,
130 positions=None, numbers=None,
131 tags=None, momenta=None, masses=None,
132 magmoms=None, charges=None,
133 scaled_positions=None,
134 cell=None, pbc=None, celldisp=None,
135 constraint=None,
136 calculator=None,
137 info=None,
138 velocities=None):
140 self._cellobj = Cell.new()
141 self._pbc = np.zeros(3, bool)
143 atoms = None
145 if hasattr(symbols, 'get_positions'):
146 atoms = symbols
147 symbols = None
148 elif (isinstance(symbols, (list, tuple)) and
149 len(symbols) > 0 and isinstance(symbols[0], Atom)):
150 # Get data from a list or tuple of Atom objects:
151 data = [[atom.get_raw(name) for atom in symbols]
152 for name in
153 ['position', 'number', 'tag', 'momentum',
154 'mass', 'magmom', 'charge']]
155 atoms = self.__class__(None, *data)
156 symbols = None
158 if atoms is not None:
159 # Get data from another Atoms object:
160 if scaled_positions is not None:
161 raise NotImplementedError
162 if symbols is None and numbers is None:
163 numbers = atoms.get_atomic_numbers()
164 if positions is None:
165 positions = atoms.get_positions()
166 if tags is None and atoms.has('tags'):
167 tags = atoms.get_tags()
168 if momenta is None and atoms.has('momenta'):
169 momenta = atoms.get_momenta()
170 if magmoms is None and atoms.has('initial_magmoms'):
171 magmoms = atoms.get_initial_magnetic_moments()
172 if masses is None and atoms.has('masses'):
173 masses = atoms.get_masses()
174 if charges is None and atoms.has('initial_charges'):
175 charges = atoms.get_initial_charges()
176 if cell is None:
177 cell = atoms.get_cell()
178 if celldisp is None:
179 celldisp = atoms.get_celldisp()
180 if pbc is None:
181 pbc = atoms.get_pbc()
182 if constraint is None:
183 constraint = [c.copy() for c in atoms.constraints]
184 if calculator is None:
185 calculator = atoms.calc
186 if info is None:
187 info = copy.deepcopy(atoms.info)
189 self.arrays = {}
191 if symbols is None:
192 if numbers is None:
193 if positions is not None:
194 natoms = len(positions)
195 elif scaled_positions is not None:
196 natoms = len(scaled_positions)
197 else:
198 natoms = 0
199 numbers = np.zeros(natoms, int)
200 self.new_array('numbers', numbers, int)
201 else:
202 if numbers is not None:
203 raise TypeError(
204 'Use only one of "symbols" and "numbers".')
205 else:
206 self.new_array('numbers', symbols2numbers(symbols), int)
208 if self.numbers.ndim != 1:
209 raise ValueError('"numbers" must be 1-dimensional.')
211 if cell is None:
212 cell = np.zeros((3, 3))
213 self.set_cell(cell)
215 if celldisp is None:
216 celldisp = np.zeros(shape=(3, 1))
217 self.set_celldisp(celldisp)
219 if positions is None:
220 if scaled_positions is None:
221 positions = np.zeros((len(self.arrays['numbers']), 3))
222 else:
223 assert self.cell.rank == 3
224 positions = np.dot(scaled_positions, self.cell)
225 else:
226 if scaled_positions is not None:
227 raise TypeError(
228 'Use only one of "symbols" and "numbers".')
229 self.new_array('positions', positions, float, (3,))
231 self.set_constraint(constraint)
232 self.set_tags(default(tags, 0))
233 self.set_masses(default(masses, None))
234 self.set_initial_magnetic_moments(default(magmoms, 0.0))
235 self.set_initial_charges(default(charges, 0.0))
236 if pbc is None:
237 pbc = False
238 self.set_pbc(pbc)
239 self.set_momenta(default(momenta, (0.0, 0.0, 0.0)),
240 apply_constraint=False)
242 if velocities is not None:
243 if momenta is None:
244 self.set_velocities(velocities)
245 else:
246 raise TypeError(
247 'Use only one of "momenta" and "velocities".')
249 if info is None:
250 self.info = {}
251 else:
252 self.info = dict(info)
254 self.calc = calculator
256 @property
257 def symbols(self):
258 """Get chemical symbols as a :class:`ase.symbols.Symbols` object.
260 The object works like ``atoms.numbers`` except its values
261 are strings. It supports in-place editing."""
262 return Symbols(self.numbers)
264 @symbols.setter
265 def symbols(self, obj):
266 new_symbols = Symbols.fromsymbols(obj)
267 self.numbers[:] = new_symbols.numbers
269 @deprecated(DeprecationWarning('Please use atoms.calc = calc'))
270 def set_calculator(self, calc=None):
271 """Attach calculator object.
273 Please use the equivalent atoms.calc = calc instead of this
274 method."""
275 self.calc = calc
277 @deprecated(DeprecationWarning('Please use atoms.calc'))
278 def get_calculator(self):
279 """Get currently attached calculator object.
281 Please use the equivalent atoms.calc instead of
282 atoms.get_calculator()."""
283 return self.calc
285 @property
286 def calc(self):
287 """Calculator object."""
288 return self._calc
290 @calc.setter
291 def calc(self, calc):
292 self._calc = calc
293 if hasattr(calc, 'set_atoms'):
294 calc.set_atoms(self)
296 @calc.deleter # type: ignore
297 @deprecated(DeprecationWarning('Please use atoms.calc = None'))
298 def calc(self):
299 self._calc = None
301 @property # type: ignore
302 @deprecated('Please use atoms.cell.rank instead')
303 def number_of_lattice_vectors(self):
304 """Number of (non-zero) lattice vectors."""
305 return self.cell.rank
307 def set_constraint(self, constraint=None):
308 """Apply one or more constrains.
310 The *constraint* argument must be one constraint object or a
311 list of constraint objects."""
312 if constraint is None:
313 self._constraints = []
314 else:
315 if isinstance(constraint, list):
316 self._constraints = constraint
317 elif isinstance(constraint, tuple):
318 self._constraints = list(constraint)
319 else:
320 self._constraints = [constraint]
322 def _get_constraints(self):
323 return self._constraints
325 def _del_constraints(self):
326 self._constraints = []
328 constraints = property(_get_constraints, set_constraint, _del_constraints,
329 'Constraints of the atoms.')
331 def set_cell(self, cell, scale_atoms=False, apply_constraint=True):
332 """Set unit cell vectors.
334 Parameters:
336 cell: 3x3 matrix or length 3 or 6 vector
337 Unit cell. A 3x3 matrix (the three unit cell vectors) or
338 just three numbers for an orthorhombic cell. Another option is
339 6 numbers, which describes unit cell with lengths of unit cell
340 vectors and with angles between them (in degrees), in following
341 order: [len(a), len(b), len(c), angle(b,c), angle(a,c),
342 angle(a,b)]. First vector will lie in x-direction, second in
343 xy-plane, and the third one in z-positive subspace.
344 scale_atoms: bool
345 Fix atomic positions or move atoms with the unit cell?
346 Default behavior is to *not* move the atoms (scale_atoms=False).
347 apply_constraint: bool
348 Whether to apply constraints to the given cell.
350 Examples:
352 Two equivalent ways to define an orthorhombic cell:
354 >>> atoms = Atoms('He')
355 >>> a, b, c = 7, 7.5, 8
356 >>> atoms.set_cell([a, b, c])
357 >>> atoms.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)])
359 FCC unit cell:
361 >>> atoms.set_cell([(0, b, b), (b, 0, b), (b, b, 0)])
363 Hexagonal unit cell:
365 >>> atoms.set_cell([a, a, c, 90, 90, 120])
367 Rhombohedral unit cell:
369 >>> alpha = 77
370 >>> atoms.set_cell([a, a, a, alpha, alpha, alpha])
371 """
373 # Override pbcs if and only if given a Cell object:
374 cell = Cell.new(cell)
376 # XXX not working well during initialize due to missing _constraints
377 if apply_constraint and hasattr(self, '_constraints'):
378 for constraint in self.constraints:
379 if hasattr(constraint, 'adjust_cell'):
380 constraint.adjust_cell(self, cell)
382 if scale_atoms:
383 M = np.linalg.solve(self.cell.complete(), cell.complete())
384 self.positions[:] = np.dot(self.positions, M)
386 self.cell[:] = cell
388 def set_celldisp(self, celldisp):
389 """Set the unit cell displacement vectors."""
390 celldisp = np.array(celldisp, float)
391 self._celldisp = celldisp
393 def get_celldisp(self):
394 """Get the unit cell displacement vectors."""
395 return self._celldisp.copy()
397 def get_cell(self, complete=False):
398 """Get the three unit cell vectors as a `class`:ase.cell.Cell` object.
400 The Cell object resembles a 3x3 ndarray, and cell[i, j]
401 is the jth Cartesian coordinate of the ith cell vector."""
402 if complete:
403 cell = self.cell.complete()
404 else:
405 cell = self.cell.copy()
407 return cell
409 @deprecated('Please use atoms.cell.cellpar() instead')
410 def get_cell_lengths_and_angles(self):
411 """Get unit cell parameters. Sequence of 6 numbers.
413 First three are unit cell vector lengths and second three
414 are angles between them::
416 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
418 in degrees.
419 """
420 return self.cell.cellpar()
422 @deprecated('Please use atoms.cell.reciprocal()')
423 def get_reciprocal_cell(self):
424 """Get the three reciprocal lattice vectors as a 3x3 ndarray.
426 Note that the commonly used factor of 2 pi for Fourier
427 transforms is not included here."""
429 return self.cell.reciprocal()
431 @property
432 def pbc(self):
433 """Reference to pbc-flags for in-place manipulations."""
434 return self._pbc
436 @pbc.setter
437 def pbc(self, pbc):
438 self._pbc[:] = pbc
440 def set_pbc(self, pbc):
441 """Set periodic boundary condition flags."""
442 self.pbc = pbc
444 def get_pbc(self):
445 """Get periodic boundary condition flags."""
446 return self.pbc.copy()
448 def new_array(self, name, a, dtype=None, shape=None):
449 """Add new array.
451 If *shape* is not *None*, the shape of *a* will be checked."""
453 if dtype is not None:
454 a = np.array(a, dtype, order='C')
455 if len(a) == 0 and shape is not None:
456 a.shape = (-1,) + shape
457 else:
458 if not a.flags['C_CONTIGUOUS']:
459 a = np.ascontiguousarray(a)
460 else:
461 a = a.copy()
463 if name in self.arrays:
464 raise RuntimeError('Array {} already present'.format(name))
466 for b in self.arrays.values():
467 if len(a) != len(b):
468 raise ValueError('Array "%s" has wrong length: %d != %d.' %
469 (name, len(a), len(b)))
470 break
472 if shape is not None and a.shape[1:] != shape:
473 raise ValueError('Array "%s" has wrong shape %s != %s.' %
474 (name, a.shape, (a.shape[0:1] + shape)))
476 self.arrays[name] = a
478 def get_array(self, name, copy=True):
479 """Get an array.
481 Returns a copy unless the optional argument copy is false.
482 """
483 if copy:
484 return self.arrays[name].copy()
485 else:
486 return self.arrays[name]
488 def set_array(self, name, a, dtype=None, shape=None):
489 """Update array.
491 If *shape* is not *None*, the shape of *a* will be checked.
492 If *a* is *None*, then the array is deleted."""
494 b = self.arrays.get(name)
495 if b is None:
496 if a is not None:
497 self.new_array(name, a, dtype, shape)
498 else:
499 if a is None:
500 del self.arrays[name]
501 else:
502 a = np.asarray(a)
503 if a.shape != b.shape:
504 raise ValueError('Array "%s" has wrong shape %s != %s.' %
505 (name, a.shape, b.shape))
506 b[:] = a
508 def has(self, name):
509 """Check for existence of array.
511 name must be one of: 'tags', 'momenta', 'masses', 'initial_magmoms',
512 'initial_charges'."""
513 # XXX extend has to calculator properties
514 return name in self.arrays
516 def set_atomic_numbers(self, numbers):
517 """Set atomic numbers."""
518 self.set_array('numbers', numbers, int, ())
520 def get_atomic_numbers(self):
521 """Get integer array of atomic numbers."""
522 return self.arrays['numbers'].copy()
524 def get_chemical_symbols(self):
525 """Get list of chemical symbol strings.
527 Equivalent to ``list(atoms.symbols)``."""
528 return list(self.symbols)
530 def set_chemical_symbols(self, symbols):
531 """Set chemical symbols."""
532 self.set_array('numbers', symbols2numbers(symbols), int, ())
534 def get_chemical_formula(self, mode='hill', empirical=False):
535 """Get the chemical formula as a string based on the chemical symbols.
537 Parameters:
539 mode: str
540 There are four different modes available:
542 'all': The list of chemical symbols are contracted to a string,
543 e.g. ['C', 'H', 'H', 'H', 'O', 'H'] becomes 'CHHHOH'.
545 'reduce': The same as 'all' where repeated elements are contracted
546 to a single symbol and a number, e.g. 'CHHHOCHHH' is reduced to
547 'CH3OCH3'.
549 'hill': The list of chemical symbols are contracted to a string
550 following the Hill notation (alphabetical order with C and H
551 first), e.g. 'CHHHOCHHH' is reduced to 'C2H6O' and 'SOOHOHO' to
552 'H2O4S'. This is default.
554 'metal': The list of chemical symbols (alphabetical metals,
555 and alphabetical non-metals)
557 empirical, bool (optional, default=False)
558 Divide the symbol counts by their greatest common divisor to yield
559 an empirical formula. Only for mode `metal` and `hill`.
560 """
561 return self.symbols.get_chemical_formula(mode, empirical)
563 def set_tags(self, tags):
564 """Set tags for all atoms. If only one tag is supplied, it is
565 applied to all atoms."""
566 if isinstance(tags, int):
567 tags = [tags] * len(self)
568 self.set_array('tags', tags, int, ())
570 def get_tags(self):
571 """Get integer array of tags."""
572 if 'tags' in self.arrays:
573 return self.arrays['tags'].copy()
574 else:
575 return np.zeros(len(self), int)
577 def set_momenta(self, momenta, apply_constraint=True):
578 """Set momenta."""
579 if (apply_constraint and len(self.constraints) > 0 and
580 momenta is not None):
581 momenta = np.array(momenta) # modify a copy
582 for constraint in self.constraints:
583 if hasattr(constraint, 'adjust_momenta'):
584 constraint.adjust_momenta(self, momenta)
585 self.set_array('momenta', momenta, float, (3,))
587 def set_velocities(self, velocities):
588 """Set the momenta by specifying the velocities."""
589 self.set_momenta(self.get_masses()[:, np.newaxis] * velocities)
591 def get_momenta(self):
592 """Get array of momenta."""
593 if 'momenta' in self.arrays:
594 return self.arrays['momenta'].copy()
595 else:
596 return np.zeros((len(self), 3))
598 def set_masses(self, masses='defaults'):
599 """Set atomic masses in atomic mass units.
601 The array masses should contain a list of masses. In case
602 the masses argument is not given or for those elements of the
603 masses list that are None, standard values are set."""
605 if isinstance(masses, str):
606 if masses == 'defaults':
607 masses = atomic_masses[self.arrays['numbers']]
608 elif masses == 'most_common':
609 masses = atomic_masses_common[self.arrays['numbers']]
610 elif masses is None:
611 pass
612 elif not isinstance(masses, np.ndarray):
613 masses = list(masses)
614 for i, mass in enumerate(masses):
615 if mass is None:
616 masses[i] = atomic_masses[self.numbers[i]]
617 self.set_array('masses', masses, float, ())
619 def get_masses(self):
620 """Get array of masses in atomic mass units."""
621 if 'masses' in self.arrays:
622 return self.arrays['masses'].copy()
623 else:
624 return atomic_masses[self.arrays['numbers']]
626 def set_initial_magnetic_moments(self, magmoms=None):
627 """Set the initial magnetic moments.
629 Use either one or three numbers for every atom (collinear
630 or non-collinear spins)."""
632 if magmoms is None:
633 self.set_array('initial_magmoms', None)
634 else:
635 magmoms = np.asarray(magmoms)
636 self.set_array('initial_magmoms', magmoms, float,
637 magmoms.shape[1:])
639 def get_initial_magnetic_moments(self):
640 """Get array of initial magnetic moments."""
641 if 'initial_magmoms' in self.arrays:
642 return self.arrays['initial_magmoms'].copy()
643 else:
644 return np.zeros(len(self))
646 def get_magnetic_moments(self):
647 """Get calculated local magnetic moments."""
648 if self._calc is None:
649 raise RuntimeError('Atoms object has no calculator.')
650 return self._calc.get_magnetic_moments(self)
652 def get_magnetic_moment(self):
653 """Get calculated total magnetic moment."""
654 if self._calc is None:
655 raise RuntimeError('Atoms object has no calculator.')
656 return self._calc.get_magnetic_moment(self)
658 def set_initial_charges(self, charges=None):
659 """Set the initial charges."""
661 if charges is None:
662 self.set_array('initial_charges', None)
663 else:
664 self.set_array('initial_charges', charges, float, ())
666 def get_initial_charges(self):
667 """Get array of initial charges."""
668 if 'initial_charges' in self.arrays:
669 return self.arrays['initial_charges'].copy()
670 else:
671 return np.zeros(len(self))
673 def get_charges(self):
674 """Get calculated charges."""
675 if self._calc is None:
676 raise RuntimeError('Atoms object has no calculator.')
677 try:
678 return self._calc.get_charges(self)
679 except AttributeError:
680 from ase.calculators.calculator import PropertyNotImplementedError
681 raise PropertyNotImplementedError
683 def set_positions(self, newpositions, apply_constraint=True):
684 """Set positions, honoring any constraints. To ignore constraints,
685 use *apply_constraint=False*."""
686 if self.constraints and apply_constraint:
687 newpositions = np.array(newpositions, float)
688 for constraint in self.constraints:
689 constraint.adjust_positions(self, newpositions)
691 self.set_array('positions', newpositions, shape=(3,))
693 def get_positions(self, wrap=False, **wrap_kw):
694 """Get array of positions.
696 Parameters:
698 wrap: bool
699 wrap atoms back to the cell before returning positions
700 wrap_kw: (keyword=value) pairs
701 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
702 see :func:`ase.geometry.wrap_positions`
703 """
704 from ase.geometry import wrap_positions
705 if wrap:
706 if 'pbc' not in wrap_kw:
707 wrap_kw['pbc'] = self.pbc
708 return wrap_positions(self.positions, self.cell, **wrap_kw)
709 else:
710 return self.arrays['positions'].copy()
712 def get_potential_energy(self, force_consistent=False,
713 apply_constraint=True):
714 """Calculate potential energy.
716 Ask the attached calculator to calculate the potential energy and
717 apply constraints. Use *apply_constraint=False* to get the raw
718 forces.
720 When supported by the calculator, either the energy extrapolated
721 to zero Kelvin or the energy consistent with the forces (the free
722 energy) can be returned.
723 """
724 if self._calc is None:
725 raise RuntimeError('Atoms object has no calculator.')
726 if force_consistent:
727 energy = self._calc.get_potential_energy(
728 self, force_consistent=force_consistent)
729 else:
730 energy = self._calc.get_potential_energy(self)
731 if apply_constraint:
732 for constraint in self.constraints:
733 if hasattr(constraint, 'adjust_potential_energy'):
734 energy += constraint.adjust_potential_energy(self)
735 return energy
737 def get_properties(self, properties):
738 """This method is experimental; currently for internal use."""
739 # XXX Something about constraints.
740 if self._calc is None:
741 raise RuntimeError('Atoms object has no calculator.')
742 return self._calc.calculate_properties(self, properties)
744 def get_potential_energies(self):
745 """Calculate the potential energies of all the atoms.
747 Only available with calculators supporting per-atom energies
748 (e.g. classical potentials).
749 """
750 if self._calc is None:
751 raise RuntimeError('Atoms object has no calculator.')
752 return self._calc.get_potential_energies(self)
754 def get_kinetic_energy(self):
755 """Get the kinetic energy."""
756 momenta = self.arrays.get('momenta')
757 if momenta is None:
758 return 0.0
759 return 0.5 * np.vdot(momenta, self.get_velocities())
761 def get_velocities(self):
762 """Get array of velocities."""
763 momenta = self.get_momenta()
764 masses = self.get_masses()
765 return momenta / masses[:, np.newaxis]
767 def get_total_energy(self):
768 """Get the total energy - potential plus kinetic energy."""
769 return self.get_potential_energy() + self.get_kinetic_energy()
771 def get_forces(self, apply_constraint=True, md=False):
772 """Calculate atomic forces.
774 Ask the attached calculator to calculate the forces and apply
775 constraints. Use *apply_constraint=False* to get the raw
776 forces.
778 For molecular dynamics (md=True) we don't apply the constraint
779 to the forces but to the momenta. When holonomic constraints for
780 rigid linear triatomic molecules are present, ask the constraints
781 to redistribute the forces within each triple defined in the
782 constraints (required for molecular dynamics with this type of
783 constraints)."""
785 if self._calc is None:
786 raise RuntimeError('Atoms object has no calculator.')
787 forces = self._calc.get_forces(self)
789 if apply_constraint:
790 # We need a special md flag here because for MD we want
791 # to skip real constraints but include special "constraints"
792 # Like Hookean.
793 for constraint in self.constraints:
794 if md and hasattr(constraint, 'redistribute_forces_md'):
795 constraint.redistribute_forces_md(self, forces)
796 if not md or hasattr(constraint, 'adjust_potential_energy'):
797 constraint.adjust_forces(self, forces)
798 return forces
800 # Informs calculators (e.g. Asap) that ideal gas contribution is added here.
801 _ase_handles_dynamic_stress = True
803 def get_stress(self, voigt=True, apply_constraint=True,
804 include_ideal_gas=False):
805 """Calculate stress tensor.
807 Returns an array of the six independent components of the
808 symmetric stress tensor, in the traditional Voigt order
809 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is Voigt
810 order.
812 The ideal gas contribution to the stresses is added if the
813 atoms have momenta and ``include_ideal_gas`` is set to True.
814 """
816 if self._calc is None:
817 raise RuntimeError('Atoms object has no calculator.')
819 stress = self._calc.get_stress(self)
820 shape = stress.shape
822 if shape == (3, 3):
823 # Convert to the Voigt form before possibly applying
824 # constraints and adding the dynamic part of the stress
825 # (the "ideal gas contribution").
826 stress = full_3x3_to_voigt_6_stress(stress)
827 else:
828 assert shape == (6,)
830 if apply_constraint:
831 for constraint in self.constraints:
832 if hasattr(constraint, 'adjust_stress'):
833 constraint.adjust_stress(self, stress)
835 # Add ideal gas contribution, if applicable
836 if include_ideal_gas and self.has('momenta'):
837 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
838 p = self.get_momenta()
839 masses = self.get_masses()
840 invmass = 1.0 / masses
841 invvol = 1.0 / self.get_volume()
842 for alpha in range(3):
843 for beta in range(alpha, 3):
844 stress[stresscomp[alpha, beta]] -= (
845 p[:, alpha] * p[:, beta] * invmass).sum() * invvol
847 if voigt:
848 return stress
849 else:
850 return voigt_6_to_full_3x3_stress(stress)
852 def get_stresses(self, include_ideal_gas=False, voigt=True):
853 """Calculate the stress-tensor of all the atoms.
855 Only available with calculators supporting per-atom energies and
856 stresses (e.g. classical potentials). Even for such calculators
857 there is a certain arbitrariness in defining per-atom stresses.
859 The ideal gas contribution to the stresses is added if the
860 atoms have momenta and ``include_ideal_gas`` is set to True.
861 """
862 if self._calc is None:
863 raise RuntimeError('Atoms object has no calculator.')
864 stresses = self._calc.get_stresses(self)
866 # make sure `stresses` are in voigt form
867 if np.shape(stresses)[1:] == (3, 3):
868 stresses_voigt = [full_3x3_to_voigt_6_stress(s) for s in stresses]
869 stresses = np.array(stresses_voigt)
871 # REMARK: The ideal gas contribution is intensive, i.e., the volume
872 # is divided out. We currently don't check if `stresses` are intensive
873 # as well, i.e., if `a.get_stresses.sum(axis=0) == a.get_stress()`.
874 # It might be good to check this here, but adds computational overhead.
876 if include_ideal_gas and self.has('momenta'):
877 stresscomp = np.array([[0, 5, 4], [5, 1, 3], [4, 3, 2]])
878 if hasattr(self._calc, 'get_atomic_volumes'):
879 invvol = 1.0 / self._calc.get_atomic_volumes()
880 else:
881 invvol = self.get_global_number_of_atoms() / self.get_volume()
882 p = self.get_momenta()
883 invmass = 1.0 / self.get_masses()
884 for alpha in range(3):
885 for beta in range(alpha, 3):
886 stresses[:, stresscomp[alpha, beta]] -= (
887 p[:, alpha] * p[:, beta] * invmass * invvol)
888 if voigt:
889 return stresses
890 else:
891 stresses_3x3 = [voigt_6_to_full_3x3_stress(s) for s in stresses]
892 return np.array(stresses_3x3)
894 def get_dipole_moment(self):
895 """Calculate the electric dipole moment for the atoms object.
897 Only available for calculators which has a get_dipole_moment()
898 method."""
900 if self._calc is None:
901 raise RuntimeError('Atoms object has no calculator.')
902 return self._calc.get_dipole_moment(self)
904 def copy(self):
905 """Return a copy."""
906 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
907 celldisp=self._celldisp.copy())
909 atoms.arrays = {}
910 for name, a in self.arrays.items():
911 atoms.arrays[name] = a.copy()
912 atoms.constraints = copy.deepcopy(self.constraints)
913 return atoms
915 def todict(self):
916 """For basic JSON (non-database) support."""
917 d = dict(self.arrays)
918 d['cell'] = np.asarray(self.cell)
919 d['pbc'] = self.pbc
920 if self._celldisp.any():
921 d['celldisp'] = self._celldisp
922 if self.constraints:
923 d['constraints'] = self.constraints
924 if self.info:
925 d['info'] = self.info
926 # Calculator... trouble.
927 return d
929 @classmethod
930 def fromdict(cls, dct):
931 """Rebuild atoms object from dictionary representation (todict)."""
932 dct = dct.copy()
933 kw = {}
934 for name in ['numbers', 'positions', 'cell', 'pbc']:
935 kw[name] = dct.pop(name)
937 constraints = dct.pop('constraints', None)
938 if constraints:
939 from ase.constraints import dict2constraint
940 constraints = [dict2constraint(d) for d in constraints]
942 info = dct.pop('info', None)
944 atoms = cls(constraint=constraints,
945 celldisp=dct.pop('celldisp', None),
946 info=info, **kw)
947 natoms = len(atoms)
949 # Some arrays are named differently from the atoms __init__ keywords.
950 # Also, there may be custom arrays. Hence we set them directly:
951 for name, arr in dct.items():
952 assert len(arr) == natoms, name
953 assert isinstance(arr, np.ndarray)
954 atoms.arrays[name] = arr
955 return atoms
957 def __len__(self):
958 return len(self.arrays['positions'])
960 def get_number_of_atoms(self):
961 """Deprecated, please do not use.
963 You probably want len(atoms). Or if your atoms are distributed,
964 use (and see) get_global_number_of_atoms()."""
965 import warnings
966 warnings.warn('Use get_global_number_of_atoms() instead',
967 np.VisibleDeprecationWarning)
968 return len(self)
970 def get_global_number_of_atoms(self):
971 """Returns the global number of atoms in a distributed-atoms parallel
972 simulation.
974 DO NOT USE UNLESS YOU KNOW WHAT YOU ARE DOING!
976 Equivalent to len(atoms) in the standard ASE Atoms class. You should
977 normally use len(atoms) instead. This function's only purpose is to
978 make compatibility between ASE and Asap easier to maintain by having a
979 few places in ASE use this function instead. It is typically only
980 when counting the global number of degrees of freedom or in similar
981 situations.
982 """
983 return len(self)
985 def __repr__(self):
986 tokens = []
988 N = len(self)
989 if N <= 60:
990 symbols = self.get_chemical_formula('reduce')
991 else:
992 symbols = self.get_chemical_formula('hill')
993 tokens.append("symbols='{0}'".format(symbols))
995 if self.pbc.any() and not self.pbc.all():
996 tokens.append('pbc={0}'.format(self.pbc.tolist()))
997 else:
998 tokens.append('pbc={0}'.format(self.pbc[0]))
1000 cell = self.cell
1001 if cell:
1002 if cell.orthorhombic:
1003 cell = cell.lengths().tolist()
1004 else:
1005 cell = cell.tolist()
1006 tokens.append('cell={0}'.format(cell))
1008 for name in sorted(self.arrays):
1009 if name in ['numbers', 'positions']:
1010 continue
1011 tokens.append('{0}=...'.format(name))
1013 if self.constraints:
1014 if len(self.constraints) == 1:
1015 constraint = self.constraints[0]
1016 else:
1017 constraint = self.constraints
1018 tokens.append('constraint={0}'.format(repr(constraint)))
1020 if self._calc is not None:
1021 tokens.append('calculator={0}(...)'
1022 .format(self._calc.__class__.__name__))
1024 return '{0}({1})'.format(self.__class__.__name__, ', '.join(tokens))
1026 def __add__(self, other):
1027 atoms = self.copy()
1028 atoms += other
1029 return atoms
1031 def extend(self, other):
1032 """Extend atoms object by appending atoms from *other*."""
1033 if isinstance(other, Atom):
1034 other = self.__class__([other])
1036 n1 = len(self)
1037 n2 = len(other)
1039 for name, a1 in self.arrays.items():
1040 a = np.zeros((n1 + n2,) + a1.shape[1:], a1.dtype)
1041 a[:n1] = a1
1042 if name == 'masses':
1043 a2 = other.get_masses()
1044 else:
1045 a2 = other.arrays.get(name)
1046 if a2 is not None:
1047 a[n1:] = a2
1048 self.arrays[name] = a
1050 for name, a2 in other.arrays.items():
1051 if name in self.arrays:
1052 continue
1053 a = np.empty((n1 + n2,) + a2.shape[1:], a2.dtype)
1054 a[n1:] = a2
1055 if name == 'masses':
1056 a[:n1] = self.get_masses()[:n1]
1057 else:
1058 a[:n1] = 0
1060 self.set_array(name, a)
1062 def __iadd__(self, other):
1063 self.extend(other)
1064 return self
1066 def append(self, atom):
1067 """Append atom to end."""
1068 self.extend(self.__class__([atom]))
1070 def __iter__(self):
1071 for i in range(len(self)):
1072 yield self[i]
1074 def __getitem__(self, i):
1075 """Return a subset of the atoms.
1077 i -- scalar integer, list of integers, or slice object
1078 describing which atoms to return.
1080 If i is a scalar, return an Atom object. If i is a list or a
1081 slice, return an Atoms object with the same cell, pbc, and
1082 other associated info as the original Atoms object. The
1083 indices of the constraints will be shuffled so that they match
1084 the indexing in the subset returned.
1086 """
1088 if isinstance(i, numbers.Integral):
1089 natoms = len(self)
1090 if i < -natoms or i >= natoms:
1091 raise IndexError('Index out of range.')
1093 return Atom(atoms=self, index=i)
1094 elif not isinstance(i, slice):
1095 i = np.array(i)
1096 if len(i) == 0:
1097 i = np.array([], dtype=int)
1098 # if i is a mask
1099 if i.dtype == bool:
1100 if len(i) != len(self):
1101 raise IndexError('Length of mask {} must equal '
1102 'number of atoms {}'
1103 .format(len(i), len(self)))
1104 i = np.arange(len(self))[i]
1106 import copy
1108 conadd = []
1109 # Constraints need to be deepcopied, but only the relevant ones.
1110 for con in copy.deepcopy(self.constraints):
1111 try:
1112 con.index_shuffle(self, i)
1113 except (IndexError, NotImplementedError):
1114 pass
1115 else:
1116 conadd.append(con)
1118 atoms = self.__class__(cell=self.cell, pbc=self.pbc, info=self.info,
1119 # should be communicated to the slice as well
1120 celldisp=self._celldisp)
1121 # TODO: Do we need to shuffle indices in adsorbate_info too?
1123 atoms.arrays = {}
1124 for name, a in self.arrays.items():
1125 atoms.arrays[name] = a[i].copy()
1127 atoms.constraints = conadd
1128 return atoms
1130 def __delitem__(self, i):
1131 from ase.constraints import FixAtoms
1132 for c in self._constraints:
1133 if not isinstance(c, FixAtoms):
1134 raise RuntimeError('Remove constraint using set_constraint() '
1135 'before deleting atoms.')
1137 if isinstance(i, list) and len(i) > 0:
1138 # Make sure a list of booleans will work correctly and not be
1139 # interpreted at 0 and 1 indices.
1140 i = np.array(i)
1142 if len(self._constraints) > 0:
1143 n = len(self)
1144 i = np.arange(n)[i]
1145 if isinstance(i, int):
1146 i = [i]
1147 constraints = []
1148 for c in self._constraints:
1149 c = c.delete_atoms(i, n)
1150 if c is not None:
1151 constraints.append(c)
1152 self.constraints = constraints
1154 mask = np.ones(len(self), bool)
1155 mask[i] = False
1156 for name, a in self.arrays.items():
1157 self.arrays[name] = a[mask]
1159 def pop(self, i=-1):
1160 """Remove and return atom at index *i* (default last)."""
1161 atom = self[i]
1162 atom.cut_reference_to_atoms()
1163 del self[i]
1164 return atom
1166 def __imul__(self, m):
1167 """In-place repeat of atoms."""
1168 if isinstance(m, int):
1169 m = (m, m, m)
1171 for x, vec in zip(m, self.cell):
1172 if x != 1 and not vec.any():
1173 raise ValueError('Cannot repeat along undefined lattice '
1174 'vector')
1176 M = np.product(m)
1177 n = len(self)
1179 for name, a in self.arrays.items():
1180 self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1))
1182 positions = self.arrays['positions']
1183 i0 = 0
1184 for m0 in range(m[0]):
1185 for m1 in range(m[1]):
1186 for m2 in range(m[2]):
1187 i1 = i0 + n
1188 positions[i0:i1] += np.dot((m0, m1, m2), self.cell)
1189 i0 = i1
1191 if self.constraints is not None:
1192 self.constraints = [c.repeat(m, n) for c in self.constraints]
1194 self.cell = np.array([m[c] * self.cell[c] for c in range(3)])
1196 return self
1198 def repeat(self, rep):
1199 """Create new repeated atoms object.
1201 The *rep* argument should be a sequence of three positive
1202 integers like *(2,3,1)* or a single integer (*r*) equivalent
1203 to *(r,r,r)*."""
1205 atoms = self.copy()
1206 atoms *= rep
1207 return atoms
1209 def __mul__(self, rep):
1210 return self.repeat(rep)
1212 def translate(self, displacement):
1213 """Translate atomic positions.
1215 The displacement argument can be a float an xyz vector or an
1216 nx3 array (where n is the number of atoms)."""
1218 self.arrays['positions'] += np.array(displacement)
1220 def center(self, vacuum=None, axis=(0, 1, 2), about=None):
1221 """Center atoms in unit cell.
1223 Centers the atoms in the unit cell, so there is the same
1224 amount of vacuum on all sides.
1226 vacuum: float (default: None)
1227 If specified adjust the amount of vacuum when centering.
1228 If vacuum=10.0 there will thus be 10 Angstrom of vacuum
1229 on each side.
1230 axis: int or sequence of ints
1231 Axis or axes to act on. Default: Act on all axes.
1232 about: float or array (default: None)
1233 If specified, center the atoms about <about>.
1234 I.e., about=(0., 0., 0.) (or just "about=0.", interpreted
1235 identically), to center about the origin.
1236 """
1238 # Find the orientations of the faces of the unit cell
1239 cell = self.cell.complete()
1240 dirs = np.zeros_like(cell)
1242 lengths = cell.lengths()
1243 for i in range(3):
1244 dirs[i] = np.cross(cell[i - 1], cell[i - 2])
1245 dirs[i] /= np.linalg.norm(dirs[i])
1246 if dirs[i] @ cell[i] < 0.0:
1247 dirs[i] *= -1
1249 if isinstance(axis, int):
1250 axes = (axis,)
1251 else:
1252 axes = axis
1254 # Now, decide how much each basis vector should be made longer
1255 pos = self.positions
1256 longer = np.zeros(3)
1257 shift = np.zeros(3)
1258 for i in axes:
1259 if len(pos):
1260 scalarprod = pos @ dirs[i]
1261 p0 = scalarprod.min()
1262 p1 = scalarprod.max()
1263 else:
1264 p0 = 0
1265 p1 = 0
1266 height = cell[i] @ dirs[i]
1267 if vacuum is not None:
1268 lng = (p1 - p0 + 2 * vacuum) - height
1269 else:
1270 lng = 0.0 # Do not change unit cell size!
1271 top = lng + height - p1
1272 shf = 0.5 * (top - p0)
1273 cosphi = cell[i] @ dirs[i] / lengths[i]
1274 longer[i] = lng / cosphi
1275 shift[i] = shf / cosphi
1277 # Now, do it!
1278 translation = np.zeros(3)
1279 for i in axes:
1280 nowlen = lengths[i]
1281 if vacuum is not None:
1282 self.cell[i] = cell[i] * (1 + longer[i] / nowlen)
1283 translation += shift[i] * cell[i] / nowlen
1285 # We calculated translations using the completed cell,
1286 # so directions without cell vectors will have been centered
1287 # along a "fake" vector of length 1.
1288 # Therefore, we adjust by -0.5:
1289 if not any(self.cell[i]):
1290 translation[i] -= 0.5
1292 # Optionally, translate to center about a point in space.
1293 if about is not None:
1294 for vector in self.cell:
1295 translation -= vector / 2.0
1296 translation += about
1298 self.positions += translation
1300 def get_center_of_mass(self, scaled=False):
1301 """Get the center of mass.
1303 If scaled=True the center of mass in scaled coordinates
1304 is returned."""
1305 masses = self.get_masses()
1306 com = masses @ self.positions / masses.sum()
1307 if scaled:
1308 return self.cell.scaled_positions(com)
1309 else:
1310 return com
1312 def set_center_of_mass(self, com, scaled=False):
1313 """Set the center of mass.
1315 If scaled=True the center of mass is expected in scaled coordinates.
1316 Constraints are considered for scaled=False.
1317 """
1318 old_com = self.get_center_of_mass(scaled=scaled)
1319 difference = com - old_com
1320 if scaled:
1321 self.set_scaled_positions(self.get_scaled_positions() + difference)
1322 else:
1323 self.set_positions(self.get_positions() + difference)
1325 def get_moments_of_inertia(self, vectors=False):
1326 """Get the moments of inertia along the principal axes.
1328 The three principal moments of inertia are computed from the
1329 eigenvalues of the symmetric inertial tensor. Periodic boundary
1330 conditions are ignored. Units of the moments of inertia are
1331 amu*angstrom**2.
1332 """
1333 com = self.get_center_of_mass()
1334 positions = self.get_positions()
1335 positions -= com # translate center of mass to origin
1336 masses = self.get_masses()
1338 # Initialize elements of the inertial tensor
1339 I11 = I22 = I33 = I12 = I13 = I23 = 0.0
1340 for i in range(len(self)):
1341 x, y, z = positions[i]
1342 m = masses[i]
1344 I11 += m * (y ** 2 + z ** 2)
1345 I22 += m * (x ** 2 + z ** 2)
1346 I33 += m * (x ** 2 + y ** 2)
1347 I12 += -m * x * y
1348 I13 += -m * x * z
1349 I23 += -m * y * z
1351 Itensor = np.array([[I11, I12, I13],
1352 [I12, I22, I23],
1353 [I13, I23, I33]])
1355 evals, evecs = np.linalg.eigh(Itensor)
1356 if vectors:
1357 return evals, evecs.transpose()
1358 else:
1359 return evals
1361 def get_angular_momentum(self):
1362 """Get total angular momentum with respect to the center of mass."""
1363 com = self.get_center_of_mass()
1364 positions = self.get_positions()
1365 positions -= com # translate center of mass to origin
1366 return np.cross(positions, self.get_momenta()).sum(0)
1368 def rotate(self, a, v, center=(0, 0, 0), rotate_cell=False):
1369 """Rotate atoms based on a vector and an angle, or two vectors.
1371 Parameters:
1373 a = None:
1374 Angle that the atoms is rotated around the vector 'v'. 'a'
1375 can also be a vector and then 'a' is rotated
1376 into 'v'.
1378 v:
1379 Vector to rotate the atoms around. Vectors can be given as
1380 strings: 'x', '-x', 'y', ... .
1382 center = (0, 0, 0):
1383 The center is kept fixed under the rotation. Use 'COM' to fix
1384 the center of mass, 'COP' to fix the center of positions or
1385 'COU' to fix the center of cell.
1387 rotate_cell = False:
1388 If true the cell is also rotated.
1390 Examples:
1392 Rotate 90 degrees around the z-axis, so that the x-axis is
1393 rotated into the y-axis:
1395 >>> atoms = Atoms()
1396 >>> atoms.rotate(90, 'z')
1397 >>> atoms.rotate(90, (0, 0, 1))
1398 >>> atoms.rotate(-90, '-z')
1399 >>> atoms.rotate('x', 'y')
1400 >>> atoms.rotate((1, 0, 0), (0, 1, 0))
1401 """
1403 if not isinstance(a, numbers.Real):
1404 a, v = v, a
1406 norm = np.linalg.norm
1407 v = string2vector(v)
1409 normv = norm(v)
1411 if normv == 0.0:
1412 raise ZeroDivisionError('Cannot rotate: norm(v) == 0')
1414 if isinstance(a, numbers.Real):
1415 a *= pi / 180
1416 v /= normv
1417 c = cos(a)
1418 s = sin(a)
1419 else:
1420 v2 = string2vector(a)
1421 v /= normv
1422 normv2 = np.linalg.norm(v2)
1423 if normv2 == 0:
1424 raise ZeroDivisionError('Cannot rotate: norm(a) == 0')
1425 v2 /= norm(v2)
1426 c = np.dot(v, v2)
1427 v = np.cross(v, v2)
1428 s = norm(v)
1429 # In case *v* and *a* are parallel, np.cross(v, v2) vanish
1430 # and can't be used as a rotation axis. However, in this
1431 # case any rotation axis perpendicular to v2 will do.
1432 eps = 1e-7
1433 if s < eps:
1434 v = np.cross((0, 0, 1), v2)
1435 if norm(v) < eps:
1436 v = np.cross((1, 0, 0), v2)
1437 assert norm(v) >= eps
1438 elif s > 0:
1439 v /= s
1441 center = self._centering_as_array(center)
1443 p = self.arrays['positions'] - center
1444 self.arrays['positions'][:] = (c * p -
1445 np.cross(p, s * v) +
1446 np.outer(np.dot(p, v), (1.0 - c) * v) +
1447 center)
1448 if rotate_cell:
1449 rotcell = self.get_cell()
1450 rotcell[:] = (c * rotcell -
1451 np.cross(rotcell, s * v) +
1452 np.outer(np.dot(rotcell, v), (1.0 - c) * v))
1453 self.set_cell(rotcell)
1455 def _centering_as_array(self, center):
1456 if isinstance(center, str):
1457 if center.lower() == 'com':
1458 center = self.get_center_of_mass()
1459 elif center.lower() == 'cop':
1460 center = self.get_positions().mean(axis=0)
1461 elif center.lower() == 'cou':
1462 center = self.get_cell().sum(axis=0) / 2
1463 else:
1464 raise ValueError('Cannot interpret center')
1465 else:
1466 center = np.array(center, float)
1467 return center
1469 def euler_rotate(self, phi=0.0, theta=0.0, psi=0.0, center=(0, 0, 0)):
1470 """Rotate atoms via Euler angles (in degrees).
1472 See e.g http://mathworld.wolfram.com/EulerAngles.html for explanation.
1474 Parameters:
1476 center :
1477 The point to rotate about. A sequence of length 3 with the
1478 coordinates, or 'COM' to select the center of mass, 'COP' to
1479 select center of positions or 'COU' to select center of cell.
1480 phi :
1481 The 1st rotation angle around the z axis.
1482 theta :
1483 Rotation around the x axis.
1484 psi :
1485 2nd rotation around the z axis.
1487 """
1488 center = self._centering_as_array(center)
1490 phi *= pi / 180
1491 theta *= pi / 180
1492 psi *= pi / 180
1494 # First move the molecule to the origin In contrast to MATLAB,
1495 # numpy broadcasts the smaller array to the larger row-wise,
1496 # so there is no need to play with the Kronecker product.
1497 rcoords = self.positions - center
1498 # First Euler rotation about z in matrix form
1499 D = np.array(((cos(phi), sin(phi), 0.),
1500 (-sin(phi), cos(phi), 0.),
1501 (0., 0., 1.)))
1502 # Second Euler rotation about x:
1503 C = np.array(((1., 0., 0.),
1504 (0., cos(theta), sin(theta)),
1505 (0., -sin(theta), cos(theta))))
1506 # Third Euler rotation, 2nd rotation about z:
1507 B = np.array(((cos(psi), sin(psi), 0.),
1508 (-sin(psi), cos(psi), 0.),
1509 (0., 0., 1.)))
1510 # Total Euler rotation
1511 A = np.dot(B, np.dot(C, D))
1512 # Do the rotation
1513 rcoords = np.dot(A, np.transpose(rcoords))
1514 # Move back to the rotation point
1515 self.positions = np.transpose(rcoords) + center
1517 def get_dihedral(self, a0, a1, a2, a3, mic=False):
1518 """Calculate dihedral angle.
1520 Calculate dihedral angle (in degrees) between the vectors a0->a1
1521 and a2->a3.
1523 Use mic=True to use the Minimum Image Convention and calculate the
1524 angle across periodic boundaries.
1525 """
1526 return self.get_dihedrals([[a0, a1, a2, a3]], mic=mic)[0]
1528 def get_dihedrals(self, indices, mic=False):
1529 """Calculate dihedral angles.
1531 Calculate dihedral angles (in degrees) between the list of vectors
1532 a0->a1 and a2->a3, where a0, a1, a2 and a3 are in each row of indices.
1534 Use mic=True to use the Minimum Image Convention and calculate the
1535 angles across periodic boundaries.
1536 """
1537 from ase.geometry import get_dihedrals
1539 indices = np.array(indices)
1540 assert indices.shape[1] == 4
1542 a0s = self.positions[indices[:, 0]]
1543 a1s = self.positions[indices[:, 1]]
1544 a2s = self.positions[indices[:, 2]]
1545 a3s = self.positions[indices[:, 3]]
1547 # vectors 0->1, 1->2, 2->3
1548 v0 = a1s - a0s
1549 v1 = a2s - a1s
1550 v2 = a3s - a2s
1552 cell = None
1553 pbc = None
1555 if mic:
1556 cell = self.cell
1557 pbc = self.pbc
1559 return get_dihedrals(v0, v1, v2, cell=cell, pbc=pbc)
1561 def _masked_rotate(self, center, axis, diff, mask):
1562 # do rotation of subgroup by copying it to temporary atoms object
1563 # and then rotating that
1564 #
1565 # recursive object definition might not be the most elegant thing,
1566 # more generally useful might be a rotation function with a mask?
1567 group = self.__class__()
1568 for i in range(len(self)):
1569 if mask[i]:
1570 group += self[i]
1571 group.translate(-center)
1572 group.rotate(diff * 180 / pi, axis)
1573 group.translate(center)
1574 # set positions in original atoms object
1575 j = 0
1576 for i in range(len(self)):
1577 if mask[i]:
1578 self.positions[i] = group[j].position
1579 j += 1
1581 def set_dihedral(self, a1, a2, a3, a4, angle,
1582 mask=None, indices=None):
1583 """Set the dihedral angle (degrees) between vectors a1->a2 and
1584 a3->a4 by changing the atom indexed by a4.
1586 If mask is not None, all the atoms described in mask
1587 (read: the entire subgroup) are moved. Alternatively to the mask,
1588 the indices of the atoms to be rotated can be supplied. If both
1589 *mask* and *indices* are given, *indices* overwrites *mask*.
1591 **Important**: If *mask* or *indices* is given and does not contain
1592 *a4*, *a4* will NOT be moved. In most cases you therefore want
1593 to include *a4* in *mask*/*indices*.
1595 Example: the following defines a very crude
1596 ethane-like molecule and twists one half of it by 30 degrees.
1598 >>> atoms = Atoms('HHCCHH', [[-1, 1, 0], [-1, -1, 0], [0, 0, 0],
1599 ... [1, 0, 0], [2, 1, 0], [2, -1, 0]])
1600 >>> atoms.set_dihedral(1, 2, 3, 4, 210, mask=[0, 0, 0, 1, 1, 1])
1601 """
1603 angle *= pi / 180
1605 # if not provided, set mask to the last atom in the
1606 # dihedral description
1607 if mask is None and indices is None:
1608 mask = np.zeros(len(self))
1609 mask[a4] = 1
1610 elif indices is not None:
1611 mask = [index in indices for index in range(len(self))]
1613 # compute necessary in dihedral change, from current value
1614 current = self.get_dihedral(a1, a2, a3, a4) * pi / 180
1615 diff = angle - current
1616 axis = self.positions[a3] - self.positions[a2]
1617 center = self.positions[a3]
1618 self._masked_rotate(center, axis, diff, mask)
1620 def rotate_dihedral(self, a1, a2, a3, a4, angle, mask=None, indices=None):
1621 """Rotate dihedral angle.
1623 Same usage as in :meth:`ase.Atoms.set_dihedral`: Rotate a group by a
1624 predefined dihedral angle, starting from its current configuration.
1625 """
1626 start = self.get_dihedral(a1, a2, a3, a4)
1627 self.set_dihedral(a1, a2, a3, a4, angle + start, mask, indices)
1629 def get_angle(self, a1, a2, a3, mic=False):
1630 """Get angle formed by three atoms.
1632 Calculate angle in degrees between the vectors a2->a1 and
1633 a2->a3.
1635 Use mic=True to use the Minimum Image Convention and calculate the
1636 angle across periodic boundaries.
1637 """
1638 return self.get_angles([[a1, a2, a3]], mic=mic)[0]
1640 def get_angles(self, indices, mic=False):
1641 """Get angle formed by three atoms for multiple groupings.
1643 Calculate angle in degrees between vectors between atoms a2->a1
1644 and a2->a3, where a1, a2, and a3 are in each row of indices.
1646 Use mic=True to use the Minimum Image Convention and calculate
1647 the angle across periodic boundaries.
1648 """
1649 from ase.geometry import get_angles
1651 indices = np.array(indices)
1652 assert indices.shape[1] == 3
1654 a1s = self.positions[indices[:, 0]]
1655 a2s = self.positions[indices[:, 1]]
1656 a3s = self.positions[indices[:, 2]]
1658 v12 = a1s - a2s
1659 v32 = a3s - a2s
1661 cell = None
1662 pbc = None
1664 if mic:
1665 cell = self.cell
1666 pbc = self.pbc
1668 return get_angles(v12, v32, cell=cell, pbc=pbc)
1670 def set_angle(self, a1, a2=None, a3=None, angle=None, mask=None,
1671 indices=None, add=False):
1672 """Set angle (in degrees) formed by three atoms.
1674 Sets the angle between vectors *a2*->*a1* and *a2*->*a3*.
1676 If *add* is `True`, the angle will be changed by the value given.
1678 Same usage as in :meth:`ase.Atoms.set_dihedral`.
1679 If *mask* and *indices*
1680 are given, *indices* overwrites *mask*. If *mask* and *indices*
1681 are not set, only *a3* is moved."""
1683 if any(a is None for a in [a2, a3, angle]):
1684 raise ValueError('a2, a3, and angle must not be None')
1686 # If not provided, set mask to the last atom in the angle description
1687 if mask is None and indices is None:
1688 mask = np.zeros(len(self))
1689 mask[a3] = 1
1690 elif indices is not None:
1691 mask = [index in indices for index in range(len(self))]
1693 if add:
1694 diff = angle
1695 else:
1696 # Compute necessary in angle change, from current value
1697 diff = angle - self.get_angle(a1, a2, a3)
1699 diff *= pi / 180
1700 # Do rotation of subgroup by copying it to temporary atoms object and
1701 # then rotating that
1702 v10 = self.positions[a1] - self.positions[a2]
1703 v12 = self.positions[a3] - self.positions[a2]
1704 v10 /= np.linalg.norm(v10)
1705 v12 /= np.linalg.norm(v12)
1706 axis = np.cross(v10, v12)
1707 center = self.positions[a2]
1708 self._masked_rotate(center, axis, diff, mask)
1710 def rattle(self, stdev=0.001, seed=None, rng=None):
1711 """Randomly displace atoms.
1713 This method adds random displacements to the atomic positions,
1714 taking a possible constraint into account. The random numbers are
1715 drawn from a normal distribution of standard deviation stdev.
1717 For a parallel calculation, it is important to use the same
1718 seed on all processors! """
1720 if seed is not None and rng is not None:
1721 raise ValueError('Please do not provide both seed and rng.')
1723 if rng is None:
1724 if seed is None:
1725 seed = 42
1726 rng = np.random.RandomState(seed)
1727 positions = self.arrays['positions']
1728 self.set_positions(positions +
1729 rng.normal(scale=stdev, size=positions.shape))
1731 def get_distance(self, a0, a1, mic=False, vector=False):
1732 """Return distance between two atoms.
1734 Use mic=True to use the Minimum Image Convention.
1735 vector=True gives the distance vector (from a0 to a1).
1736 """
1737 return self.get_distances(a0, [a1], mic=mic, vector=vector)[0]
1739 def get_distances(self, a, indices, mic=False, vector=False):
1740 """Return distances of atom No.i with a list of atoms.
1742 Use mic=True to use the Minimum Image Convention.
1743 vector=True gives the distance vector (from a to self[indices]).
1744 """
1745 from ase.geometry import get_distances
1747 R = self.arrays['positions']
1748 p1 = [R[a]]
1749 p2 = R[indices]
1751 cell = None
1752 pbc = None
1754 if mic:
1755 cell = self.cell
1756 pbc = self.pbc
1758 D, D_len = get_distances(p1, p2, cell=cell, pbc=pbc)
1760 if vector:
1761 D.shape = (-1, 3)
1762 return D
1763 else:
1764 D_len.shape = (-1,)
1765 return D_len
1767 def get_all_distances(self, mic=False, vector=False):
1768 """Return distances of all of the atoms with all of the atoms.
1770 Use mic=True to use the Minimum Image Convention.
1771 """
1772 from ase.geometry import get_distances
1774 R = self.arrays['positions']
1776 cell = None
1777 pbc = None
1779 if mic:
1780 cell = self.cell
1781 pbc = self.pbc
1783 D, D_len = get_distances(R, cell=cell, pbc=pbc)
1785 if vector:
1786 return D
1787 else:
1788 return D_len
1790 def set_distance(self, a0, a1, distance, fix=0.5, mic=False,
1791 mask=None, indices=None, add=False, factor=False):
1792 """Set the distance between two atoms.
1794 Set the distance between atoms *a0* and *a1* to *distance*.
1795 By default, the center of the two atoms will be fixed. Use
1796 *fix=0* to fix the first atom, *fix=1* to fix the second
1797 atom and *fix=0.5* (default) to fix the center of the bond.
1799 If *mask* or *indices* are set (*mask* overwrites *indices*),
1800 only the atoms defined there are moved
1801 (see :meth:`ase.Atoms.set_dihedral`).
1803 When *add* is true, the distance is changed by the value given.
1804 In combination
1805 with *factor* True, the value given is a factor scaling the distance.
1807 It is assumed that the atoms in *mask*/*indices* move together
1808 with *a1*. If *fix=1*, only *a0* will therefore be moved."""
1809 from ase.geometry import find_mic
1811 if a0 % len(self) == a1 % len(self):
1812 raise ValueError('a0 and a1 must not be the same')
1814 if add:
1815 oldDist = self.get_distance(a0, a1, mic=mic)
1816 if factor:
1817 newDist = oldDist * distance
1818 else:
1819 newDist = oldDist + distance
1820 self.set_distance(a0, a1, newDist, fix=fix, mic=mic,
1821 mask=mask, indices=indices, add=False,
1822 factor=False)
1823 return
1825 R = self.arrays['positions']
1826 D = np.array([R[a1] - R[a0]])
1828 if mic:
1829 D, D_len = find_mic(D, self.cell, self.pbc)
1830 else:
1831 D_len = np.array([np.sqrt((D**2).sum())])
1832 x = 1.0 - distance / D_len[0]
1834 if mask is None and indices is None:
1835 indices = [a0, a1]
1836 elif mask:
1837 indices = [i for i in range(len(self)) if mask[i]]
1839 for i in indices:
1840 if i == a0:
1841 R[a0] += (x * fix) * D[0]
1842 else:
1843 R[i] -= (x * (1.0 - fix)) * D[0]
1845 def get_scaled_positions(self, wrap=True):
1846 """Get positions relative to unit cell.
1848 If wrap is True, atoms outside the unit cell will be wrapped into
1849 the cell in those directions with periodic boundary conditions
1850 so that the scaled coordinates are between zero and one.
1852 If any cell vectors are zero, the corresponding coordinates
1853 are evaluated as if the cell were completed using
1854 ``cell.complete()``. This means coordinates will be Cartesian
1855 as long as the non-zero cell vectors span a Cartesian axis or
1856 plane."""
1858 fractional = self.cell.scaled_positions(self.positions)
1860 if wrap:
1861 for i, periodic in enumerate(self.pbc):
1862 if periodic:
1863 # Yes, we need to do it twice.
1864 # See the scaled_positions.py test.
1865 fractional[:, i] %= 1.0
1866 fractional[:, i] %= 1.0
1868 return fractional
1870 def set_scaled_positions(self, scaled):
1871 """Set positions relative to unit cell."""
1872 self.positions[:] = self.cell.cartesian_positions(scaled)
1874 def wrap(self, **wrap_kw):
1875 """Wrap positions to unit cell.
1877 Parameters:
1879 wrap_kw: (keyword=value) pairs
1880 optional keywords `pbc`, `center`, `pretty_translation`, `eps`,
1881 see :func:`ase.geometry.wrap_positions`
1882 """
1884 if 'pbc' not in wrap_kw:
1885 wrap_kw['pbc'] = self.pbc
1887 self.positions[:] = self.get_positions(wrap=True, **wrap_kw)
1889 def get_temperature(self):
1890 """Get the temperature in Kelvin."""
1891 dof = len(self) * 3
1892 for constraint in self._constraints:
1893 dof -= constraint.get_removed_dof(self)
1894 ekin = self.get_kinetic_energy()
1895 return 2 * ekin / (dof * units.kB)
1897 def __eq__(self, other):
1898 """Check for identity of two atoms objects.
1900 Identity means: same positions, atomic numbers, unit cell and
1901 periodic boundary conditions."""
1902 if not isinstance(other, Atoms):
1903 return False
1904 a = self.arrays
1905 b = other.arrays
1906 return (len(self) == len(other) and
1907 (a['positions'] == b['positions']).all() and
1908 (a['numbers'] == b['numbers']).all() and
1909 (self.cell == other.cell).all() and
1910 (self.pbc == other.pbc).all())
1912 def __ne__(self, other):
1913 """Check if two atoms objects are not equal.
1915 Any differences in positions, atomic numbers, unit cell or
1916 periodic boundary condtions make atoms objects not equal.
1917 """
1918 eq = self.__eq__(other)
1919 if eq is NotImplemented:
1920 return eq
1921 else:
1922 return not eq
1924 # @deprecated('Please use atoms.cell.volume')
1925 # We kind of want to deprecate this, but the ValueError behaviour
1926 # might be desirable. Should we do this?
1927 def get_volume(self):
1928 """Get volume of unit cell."""
1929 if self.cell.rank != 3:
1930 raise ValueError(
1931 'You have {0} lattice vectors: volume not defined'
1932 .format(self.cell.rank))
1933 return self.cell.volume
1935 def _get_positions(self):
1936 """Return reference to positions-array for in-place manipulations."""
1937 return self.arrays['positions']
1939 def _set_positions(self, pos):
1940 """Set positions directly, bypassing constraints."""
1941 self.arrays['positions'][:] = pos
1943 positions = property(_get_positions, _set_positions,
1944 doc='Attribute for direct ' +
1945 'manipulation of the positions.')
1947 def _get_atomic_numbers(self):
1948 """Return reference to atomic numbers for in-place
1949 manipulations."""
1950 return self.arrays['numbers']
1952 numbers = property(_get_atomic_numbers, set_atomic_numbers,
1953 doc='Attribute for direct ' +
1954 'manipulation of the atomic numbers.')
1956 @property
1957 def cell(self):
1958 """The :class:`ase.cell.Cell` for direct manipulation."""
1959 return self._cellobj
1961 @cell.setter
1962 def cell(self, cell):
1963 cell = Cell.ascell(cell)
1964 self._cellobj[:] = cell
1966 def write(self, filename, format=None, **kwargs):
1967 """Write atoms object to a file.
1969 see ase.io.write for formats.
1970 kwargs are passed to ase.io.write.
1971 """
1972 from ase.io import write
1973 write(filename, self, format, **kwargs)
1975 def iterimages(self):
1976 yield self
1978 def edit(self):
1979 """Modify atoms interactively through ASE's GUI viewer.
1981 Conflicts leading to undesirable behaviour might arise
1982 when matplotlib has been pre-imported with certain
1983 incompatible backends and while trying to use the
1984 plot feature inside the interactive GUI. To circumvent,
1985 please set matplotlib.use('gtk') before calling this
1986 method.
1987 """
1988 from ase.gui.images import Images
1989 from ase.gui.gui import GUI
1990 images = Images([self])
1991 gui = GUI(images)
1992 gui.run()
1995def string2vector(v):
1996 if isinstance(v, str):
1997 if v[0] == '-':
1998 return -string2vector(v[1:])
1999 w = np.zeros(3)
2000 w['xyz'.index(v)] = 1.0
2001 return w
2002 return np.array(v, float)
2005def default(data, dflt):
2006 """Helper function for setting default values."""
2007 if data is None:
2008 return None
2009 elif isinstance(data, (list, tuple)):
2010 newdata = []
2011 allnone = True
2012 for x in data:
2013 if x is None:
2014 newdata.append(dflt)
2015 else:
2016 newdata.append(x)
2017 allnone = False
2018 if allnone:
2019 return None
2020 return newdata
2021 else:
2022 return data