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