Coverage for /builds/alexhroom/ase/ase/constraints.py: 87.12%
1219 statements
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-05 14:37 +0000
1"""Constraints"""
2from typing import Sequence
3from warnings import warn
5import numpy as np
7from ase import Atoms
8from ase.filters import ExpCellFilter as ExpCellFilterOld
9from ase.filters import Filter as FilterOld
10from ase.filters import StrainFilter as StrainFilterOld
11from ase.filters import UnitCellFilter as UnitCellFilterOld
12from ase.geometry import (conditional_find_mic, find_mic, get_angles,
13 get_angles_derivatives, get_dihedrals,
14 get_dihedrals_derivatives, get_distances_derivatives,
15 wrap_positions)
16from ase.spacegroup.symmetrize import (prep_symmetry, refine_symmetry,
17 symmetrize_rank1, symmetrize_rank2)
18from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
19from ase.utils import deprecated
20from ase.utils.parsemath import eval_expression
22__all__ = [
23 'FixCartesian', 'FixBondLength', 'FixedMode',
24 'FixAtoms', 'FixScaled', 'FixCom', 'FixSubsetCom', 'FixedPlane',
25 'FixConstraint', 'FixedLine', 'FixBondLengths', 'FixLinearTriatomic',
26 'FixInternals', 'Hookean', 'ExternalForce', 'MirrorForce', 'MirrorTorque',
27 'FixScaledParametricRelations', 'FixCartesianParametricRelations',
28 'FixSymmetry']
31def dict2constraint(dct):
32 if dct['name'] not in __all__:
33 raise ValueError
34 return globals()[dct['name']](**dct['kwargs'])
37def slice2enlist(s, n):
38 """Convert a slice object into a list of (new, old) tuples."""
39 if isinstance(s, slice):
40 return enumerate(range(*s.indices(n)))
41 return enumerate(s)
44def constrained_indices(atoms, only_include=None):
45 """Returns a list of indices for the atoms that are constrained
46 by a constraint that is applied. By setting only_include to a
47 specific type of constraint you can make it only look for that
48 given constraint.
49 """
50 indices = []
51 for constraint in atoms.constraints:
52 if only_include is not None:
53 if not isinstance(constraint, only_include):
54 continue
55 indices.extend(np.array(constraint.get_indices()))
56 return np.array(np.unique(indices))
59class FixConstraint:
60 """Base class for classes that fix one or more atoms in some way."""
62 def index_shuffle(self, atoms: Atoms, ind):
63 """Change the indices.
65 When the ordering of the atoms in the Atoms object changes,
66 this method can be called to shuffle the indices of the
67 constraints.
69 ind -- List or tuple of indices.
71 """
72 raise NotImplementedError
74 def repeat(self, m: int, n: int):
75 """ basic method to multiply by m, needs to know the length
76 of the underlying atoms object for the assignment of
77 multiplied constraints to work.
78 """
79 msg = ("Repeat is not compatible with your atoms' constraints."
80 ' Use atoms.set_constraint() before calling repeat to '
81 'remove your constraints.')
82 raise NotImplementedError(msg)
84 def get_removed_dof(self, atoms: Atoms):
85 """Get number of removed degrees of freedom due to constraint."""
86 raise NotImplementedError
88 def adjust_positions(self, atoms: Atoms, new):
89 """Adjust positions."""
91 def adjust_momenta(self, atoms: Atoms, momenta):
92 """Adjust momenta."""
93 # The default is in identical manner to forces.
94 # TODO: The default is however not always reasonable.
95 self.adjust_forces(atoms, momenta)
97 def adjust_forces(self, atoms: Atoms, forces):
98 """Adjust forces."""
100 def copy(self):
101 """Copy constraint."""
102 return dict2constraint(self.todict().copy())
104 def todict(self):
105 """Convert constraint to dictionary."""
108class IndexedConstraint(FixConstraint):
109 def __init__(self, indices=None, mask=None):
110 """Constrain chosen atoms.
112 Parameters
113 ----------
114 indices : sequence of int
115 Indices for those atoms that should be constrained.
116 mask : sequence of bool
117 One boolean per atom indicating if the atom should be
118 constrained or not.
119 """
121 if mask is not None:
122 if indices is not None:
123 raise ValueError('Use only one of "indices" and "mask".')
124 indices = mask
125 indices = np.atleast_1d(indices)
126 if np.ndim(indices) > 1:
127 raise ValueError('indices has wrong amount of dimensions. '
128 f'Got {np.ndim(indices)}, expected ndim <= 1')
130 if indices.dtype == bool:
131 indices = np.arange(len(indices))[indices]
132 elif len(indices) == 0:
133 indices = np.empty(0, dtype=int)
134 elif not np.issubdtype(indices.dtype, np.integer):
135 raise ValueError('Indices must be integers or boolean mask, '
136 f'not dtype={indices.dtype}')
138 if len(set(indices)) < len(indices):
139 raise ValueError(
140 'The indices array contains duplicates. '
141 'Perhaps you want to specify a mask instead, but '
142 'forgot the mask= keyword.')
144 self.index = indices
146 def index_shuffle(self, atoms, ind):
147 # See docstring of superclass
148 index = []
150 # Resolve negative indices:
151 actual_indices = set(np.arange(len(atoms))[self.index])
153 for new, old in slice2enlist(ind, len(atoms)):
154 if old in actual_indices:
155 index.append(new)
156 if len(index) == 0:
157 raise IndexError('All indices in FixAtoms not part of slice')
158 self.index = np.asarray(index, int)
159 # XXX make immutable
161 def get_indices(self):
162 return self.index.copy()
164 def repeat(self, m, n):
165 i0 = 0
166 natoms = 0
167 if isinstance(m, int):
168 m = (m, m, m)
169 index_new = []
170 for _ in range(m[2]):
171 for _ in range(m[1]):
172 for _ in range(m[0]):
173 i1 = i0 + n
174 index_new += [i + natoms for i in self.index]
175 i0 = i1
176 natoms += n
177 self.index = np.asarray(index_new, int)
178 # XXX make immutable
179 return self
181 def delete_atoms(self, indices, natoms):
182 """Removes atoms from the index array, if present.
184 Required for removing atoms with existing constraint.
185 """
187 i = np.zeros(natoms, int) - 1
188 new = np.delete(np.arange(natoms), indices)
189 i[new] = np.arange(len(new))
190 index = i[self.index]
191 self.index = index[index >= 0]
192 # XXX make immutable
193 if len(self.index) == 0:
194 return None
195 return self
198class FixAtoms(IndexedConstraint):
199 """Fix chosen atoms.
201 Examples
202 --------
203 Fix all Copper atoms:
205 >>> from ase.build import bulk
207 >>> atoms = bulk('Cu', 'fcc', a=3.6)
208 >>> mask = (atoms.symbols == 'Cu')
209 >>> c = FixAtoms(mask=mask)
210 >>> atoms.set_constraint(c)
212 Fix all atoms with z-coordinate less than 1.0 Angstrom:
214 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)
215 >>> atoms.set_constraint(c)
216 """
218 def get_removed_dof(self, atoms):
219 return 3 * len(self.index)
221 def adjust_positions(self, atoms, new):
222 new[self.index] = atoms.positions[self.index]
224 def adjust_forces(self, atoms, forces):
225 forces[self.index] = 0.0
227 def __repr__(self):
228 clsname = type(self).__name__
229 indices = ints2string(self.index)
230 return f'{clsname}(indices={indices})'
232 def todict(self):
233 return {'name': 'FixAtoms',
234 'kwargs': {'indices': self.index.tolist()}}
237class FixCom(FixConstraint):
238 """Constraint class for fixing the center of mass."""
240 index = slice(None) # all atoms
242 def get_removed_dof(self, atoms):
243 return 3
245 def adjust_positions(self, atoms, new):
246 masses = atoms.get_masses()[self.index]
247 old_cm = atoms.get_center_of_mass(indices=self.index)
248 new_cm = masses @ new[self.index] / masses.sum()
249 diff = old_cm - new_cm
250 new += diff
252 def adjust_momenta(self, atoms, momenta):
253 """Adjust momenta so that the center-of-mass velocity is zero."""
254 masses = atoms.get_masses()[self.index]
255 velocity_com = momenta[self.index].sum(axis=0) / masses.sum()
256 momenta[self.index] -= masses[:, None] * velocity_com
258 def adjust_forces(self, atoms, forces):
259 # Eqs. (3) and (7) in https://doi.org/10.1021/jp9722824
260 masses = atoms.get_masses()[self.index]
261 lmd = masses @ forces[self.index] / sum(masses**2)
262 forces[self.index] -= masses[:, None] * lmd
264 def todict(self):
265 return {'name': 'FixCom',
266 'kwargs': {}}
269class FixSubsetCom(FixCom, IndexedConstraint):
270 """Constraint class for fixing the center of mass of a subset of atoms."""
272 def __init__(self, indices):
273 super().__init__(indices=indices)
275 def todict(self):
276 return {'name': self.__class__.__name__,
277 'kwargs': {'indices': self.index.tolist()}}
280def ints2string(x, threshold=None):
281 """Convert ndarray of ints to string."""
282 if threshold is None or len(x) <= threshold:
283 return str(x.tolist())
284 return str(x[:threshold].tolist())[:-1] + ', ...]'
287class FixBondLengths(FixConstraint):
288 maxiter = 500
290 def __init__(self, pairs, tolerance=1e-13,
291 bondlengths=None, iterations=None):
292 """iterations:
293 Ignored"""
294 self.pairs = np.asarray(pairs)
295 self.tolerance = tolerance
296 self.bondlengths = bondlengths
298 def get_removed_dof(self, atoms):
299 return len(self.pairs)
301 def adjust_positions(self, atoms, new):
302 old = atoms.positions
303 masses = atoms.get_masses()
305 if self.bondlengths is None:
306 self.bondlengths = self.initialize_bond_lengths(atoms)
308 for i in range(self.maxiter):
309 converged = True
310 for j, ab in enumerate(self.pairs):
311 a = ab[0]
312 b = ab[1]
313 cd = self.bondlengths[j]
314 r0 = old[a] - old[b]
315 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
316 d1 = new[a] - new[b] - r0 + d0
317 m = 1 / (1 / masses[a] + 1 / masses[b])
318 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)
319 if abs(x) > self.tolerance:
320 new[a] += x * m / masses[a] * d0
321 new[b] -= x * m / masses[b] * d0
322 converged = False
323 if converged:
324 break
325 else:
326 raise RuntimeError('Did not converge')
328 def adjust_momenta(self, atoms, p):
329 old = atoms.positions
330 masses = atoms.get_masses()
332 if self.bondlengths is None:
333 self.bondlengths = self.initialize_bond_lengths(atoms)
335 for i in range(self.maxiter):
336 converged = True
337 for j, ab in enumerate(self.pairs):
338 a = ab[0]
339 b = ab[1]
340 cd = self.bondlengths[j]
341 d = old[a] - old[b]
342 d, _ = find_mic(d, atoms.cell, atoms.pbc)
343 dv = p[a] / masses[a] - p[b] / masses[b]
344 m = 1 / (1 / masses[a] + 1 / masses[b])
345 x = -np.dot(dv, d) / cd**2
346 if abs(x) > self.tolerance:
347 p[a] += x * m * d
348 p[b] -= x * m * d
349 converged = False
350 if converged:
351 break
352 else:
353 raise RuntimeError('Did not converge')
355 def adjust_forces(self, atoms, forces):
356 self.constraint_forces = -forces
357 self.adjust_momenta(atoms, forces)
358 self.constraint_forces += forces
360 def initialize_bond_lengths(self, atoms):
361 bondlengths = np.zeros(len(self.pairs))
363 for i, ab in enumerate(self.pairs):
364 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)
366 return bondlengths
368 def get_indices(self):
369 return np.unique(self.pairs.ravel())
371 def todict(self):
372 return {'name': 'FixBondLengths',
373 'kwargs': {'pairs': self.pairs.tolist(),
374 'tolerance': self.tolerance}}
376 def index_shuffle(self, atoms, ind):
377 """Shuffle the indices of the two atoms in this constraint"""
378 map = np.zeros(len(atoms), int)
379 map[ind] = 1
380 n = map.sum()
381 map[:] = -1
382 map[ind] = range(n)
383 pairs = map[self.pairs]
384 self.pairs = pairs[(pairs != -1).all(1)]
385 if len(self.pairs) == 0:
386 raise IndexError('Constraint not part of slice')
389def FixBondLength(a1, a2):
390 """Fix distance between atoms with indices a1 and a2."""
391 return FixBondLengths([(a1, a2)])
394class FixLinearTriatomic(FixConstraint):
395 """Holonomic constraints for rigid linear triatomic molecules."""
397 def __init__(self, triples):
398 """Apply RATTLE-type bond constraints between outer atoms n and m
399 and linear vectorial constraints to the position of central
400 atoms o to fix the geometry of linear triatomic molecules of the
401 type:
403 n--o--m
405 Parameters:
407 triples: list
408 Indices of the atoms forming the linear molecules to constrain
409 as triples. Sequence should be (n, o, m) or (m, o, n).
411 When using these constraints in molecular dynamics or structure
412 optimizations, atomic forces need to be redistributed within a
413 triple. The function redistribute_forces_optimization implements
414 the redistribution of forces for structure optimization, while
415 the function redistribute_forces_md implements the redistribution
416 for molecular dynamics.
418 References:
420 Ciccotti et al. Molecular Physics 47 (1982)
421 :doi:`10.1080/00268978200100942`
422 """
423 self.triples = np.asarray(triples)
424 if self.triples.shape[1] != 3:
425 raise ValueError('"triples" has wrong size')
426 self.bondlengths = None
428 def get_removed_dof(self, atoms):
429 return 4 * len(self.triples)
431 @property
432 def n_ind(self):
433 return self.triples[:, 0]
435 @property
436 def m_ind(self):
437 return self.triples[:, 2]
439 @property
440 def o_ind(self):
441 return self.triples[:, 1]
443 def initialize(self, atoms):
444 masses = atoms.get_masses()
445 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)
447 self.bondlengths = self.initialize_bond_lengths(atoms)
448 self.bondlengths_nm = self.bondlengths.sum(axis=1)
450 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]
451 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +
452 C1[:, 1] ** 2 * self.mass_n * self.mass_o +
453 self.mass_n * self.mass_m)
454 C2 = C1 / C2[:, None]
455 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]
456 C3 = C2 * self.mass_o[:, None] * C3[:, None]
457 C3[:, 1] *= -1
458 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T
459 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)
460 C4 = C1 / C4[:, None]
462 self.C1 = C1
463 self.C2 = C2
464 self.C3 = C3
465 self.C4 = C4
467 def adjust_positions(self, atoms, new):
468 old = atoms.positions
469 new_n, new_m, new_o = self.get_slices(new)
471 if self.bondlengths is None:
472 self.initialize(atoms)
474 r0 = old[self.n_ind] - old[self.m_ind]
475 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
476 d1 = new_n - new_m - r0 + d0
477 a = np.einsum('ij,ij->i', d0, d0)
478 b = np.einsum('ij,ij->i', d1, d0)
479 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2
480 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))
481 g = g[:, None] * self.C3
482 new_n -= g[:, 0, None] * d0
483 new_m += g[:, 1, None] * d0
484 if np.allclose(d0, r0):
485 new_o = (self.C1[:, 0, None] * new_n
486 + self.C1[:, 1, None] * new_m)
487 else:
488 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)
489 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)
490 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2
491 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)
493 self.set_slices(new_n, new_m, new_o, new)
495 def adjust_momenta(self, atoms, p):
496 old = atoms.positions
497 p_n, p_m, p_o = self.get_slices(p)
499 if self.bondlengths is None:
500 self.initialize(atoms)
502 mass_nn = self.mass_n[:, None]
503 mass_mm = self.mass_m[:, None]
504 mass_oo = self.mass_o[:, None]
506 d = old[self.n_ind] - old[self.m_ind]
507 d, _ = find_mic(d, atoms.cell, atoms.pbc)
508 dv = p_n / mass_nn - p_m / mass_mm
509 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2
510 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]
511 p_n -= k[:, 0, None] * mass_nn * d
512 p_m += k[:, 1, None] * mass_mm * d
513 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +
514 self.C1[:, 1, None] * p_m / mass_mm))
516 self.set_slices(p_n, p_m, p_o, p)
518 def adjust_forces(self, atoms, forces):
520 if self.bondlengths is None:
521 self.initialize(atoms)
523 A = self.C4 * np.diff(self.C1)
524 A[:, 0] *= -1
525 A -= 1
526 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]
527 A /= (A.sum(axis=1))[:, None]
529 self.constraint_forces = -forces
530 old = atoms.positions
532 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)
534 d = old[self.n_ind] - old[self.m_ind]
535 d, _ = find_mic(d, atoms.cell, atoms.pbc)
536 df = fr_n - fr_m
537 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2
538 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]
539 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]
540 forces[self.o_ind] = fr_o + k[:, None] * d * B
542 self.constraint_forces += forces
544 def redistribute_forces_optimization(self, forces):
545 """Redistribute forces within a triple when performing structure
546 optimizations.
548 The redistributed forces needs to be further adjusted using the
549 appropriate Lagrange multipliers as implemented in adjust_forces."""
550 forces_n, forces_m, forces_o = self.get_slices(forces)
551 C1_1 = self.C1[:, 0, None]
552 C1_2 = self.C1[:, 1, None]
553 C4_1 = self.C4[:, 0, None]
554 C4_2 = self.C4[:, 1, None]
556 fr_n = ((1 - C4_1 * C1_1) * forces_n -
557 C4_1 * (C1_2 * forces_m - forces_o))
558 fr_m = ((1 - C4_2 * C1_2) * forces_m -
559 C4_2 * (C1_1 * forces_n - forces_o))
560 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +
561 C4_1 * forces_n + C4_2 * forces_m)
563 return fr_n, fr_m, fr_o
565 def redistribute_forces_md(self, atoms, forces, rand=False):
566 """Redistribute forces within a triple when performing molecular
567 dynamics.
569 When rand=True, use the equations for random force terms, as
570 used e.g. by Langevin dynamics, otherwise apply the standard
571 equations for deterministic forces (see Ciccotti et al. Molecular
572 Physics 47 (1982))."""
573 if self.bondlengths is None:
574 self.initialize(atoms)
575 forces_n, forces_m, forces_o = self.get_slices(forces)
576 C1_1 = self.C1[:, 0, None]
577 C1_2 = self.C1[:, 1, None]
578 C2_1 = self.C2[:, 0, None]
579 C2_2 = self.C2[:, 1, None]
580 mass_nn = self.mass_n[:, None]
581 mass_mm = self.mass_m[:, None]
582 mass_oo = self.mass_o[:, None]
583 if rand:
584 mr1 = (mass_mm / mass_nn) ** 0.5
585 mr2 = (mass_oo / mass_nn) ** 0.5
586 mr3 = (mass_nn / mass_mm) ** 0.5
587 mr4 = (mass_oo / mass_mm) ** 0.5
588 else:
589 mr1 = 1.0
590 mr2 = 1.0
591 mr3 = 1.0
592 mr4 = 1.0
594 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -
595 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -
596 mr2 * mass_mm * mass_nn * forces_o))
598 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -
599 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -
600 mr4 * mass_mm * mass_nn * forces_o))
602 self.set_slices(fr_n, fr_m, 0.0, forces)
604 def get_slices(self, a):
605 a_n = a[self.n_ind]
606 a_m = a[self.m_ind]
607 a_o = a[self.o_ind]
609 return a_n, a_m, a_o
611 def set_slices(self, a_n, a_m, a_o, a):
612 a[self.n_ind] = a_n
613 a[self.m_ind] = a_m
614 a[self.o_ind] = a_o
616 def initialize_bond_lengths(self, atoms):
617 bondlengths = np.zeros((len(self.triples), 2))
619 for i in range(len(self.triples)):
620 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],
621 self.o_ind[i], mic=True)
622 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],
623 self.m_ind[i], mic=True)
625 return bondlengths
627 def get_indices(self):
628 return np.unique(self.triples.ravel())
630 def todict(self):
631 return {'name': 'FixLinearTriatomic',
632 'kwargs': {'triples': self.triples.tolist()}}
634 def index_shuffle(self, atoms, ind):
635 """Shuffle the indices of the three atoms in this constraint"""
636 map = np.zeros(len(atoms), int)
637 map[ind] = 1
638 n = map.sum()
639 map[:] = -1
640 map[ind] = range(n)
641 triples = map[self.triples]
642 self.triples = triples[(triples != -1).all(1)]
643 if len(self.triples) == 0:
644 raise IndexError('Constraint not part of slice')
647class FixedMode(FixConstraint):
648 """Constrain atoms to move along directions orthogonal to
649 a given mode only. Initialize with a mode, such as one produced by
650 ase.vibrations.Vibrations.get_mode()."""
652 def __init__(self, mode):
653 mode = np.asarray(mode)
654 self.mode = (mode / np.sqrt((mode**2).sum())).reshape(-1)
656 def get_removed_dof(self, atoms):
657 return len(atoms)
659 def adjust_positions(self, atoms, newpositions):
660 newpositions = newpositions.ravel()
661 oldpositions = atoms.positions.ravel()
662 step = newpositions - oldpositions
663 newpositions -= self.mode * np.dot(step, self.mode)
665 def adjust_forces(self, atoms, forces):
666 forces = forces.ravel()
667 forces -= self.mode * np.dot(forces, self.mode)
669 def index_shuffle(self, atoms, ind):
670 eps = 1e-12
671 mode = self.mode.reshape(-1, 3)
672 excluded = np.ones(len(mode), dtype=bool)
673 excluded[ind] = False
674 if (abs(mode[excluded]) > eps).any():
675 raise IndexError('All nonzero parts of mode not in slice')
676 self.mode = mode[ind].ravel()
678 def get_indices(self):
679 # This function will never properly work because it works on all
680 # atoms and it has no idea how to tell how many atoms it is
681 # attached to. If it is being used, surely the user knows
682 # everything is being constrained.
683 return []
685 def todict(self):
686 return {'name': 'FixedMode',
687 'kwargs': {'mode': self.mode.tolist()}}
689 def __repr__(self):
690 return f'FixedMode({self.mode.tolist()})'
693def _normalize(direction):
694 if np.shape(direction) != (3,):
695 raise ValueError("len(direction) is {len(direction)}. Has to be 3")
697 direction = np.asarray(direction) / np.linalg.norm(direction)
698 return direction
701class FixedPlane(IndexedConstraint):
702 """
703 Constraint object for fixing chosen atoms to only move in a plane.
705 The plane is defined by its normal vector *direction*
706 """
708 def __init__(self, indices, direction):
709 """Constrain chosen atoms.
711 Parameters
712 ----------
713 indices : int or list of int
714 Index or indices for atoms that should be constrained
715 direction : list of 3 int
716 Direction of the normal vector
718 Examples
719 --------
720 Fix all Copper atoms to only move in the yz-plane:
722 >>> from ase.build import bulk
723 >>> from ase.constraints import FixedPlane
725 >>> atoms = bulk('Cu', 'fcc', a=3.6)
726 >>> c = FixedPlane(
727 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
728 ... direction=[1, 0, 0],
729 ... )
730 >>> atoms.set_constraint(c)
732 or constrain a single atom with the index 0 to move in the xy-plane:
734 >>> c = FixedPlane(indices=0, direction=[0, 0, 1])
735 >>> atoms.set_constraint(c)
736 """
737 super().__init__(indices=indices)
738 self.dir = _normalize(direction)
740 def adjust_positions(self, atoms, newpositions):
741 step = newpositions[self.index] - atoms.positions[self.index]
742 newpositions[self.index] -= _projection(step, self.dir)
744 def adjust_forces(self, atoms, forces):
745 forces[self.index] -= _projection(forces[self.index], self.dir)
747 def get_removed_dof(self, atoms):
748 return len(self.index)
750 def todict(self):
751 return {
752 'name': 'FixedPlane',
753 'kwargs': {'indices': self.index.tolist(),
754 'direction': self.dir.tolist()}
755 }
757 def __repr__(self):
758 return f'FixedPlane(indices={self.index}, {self.dir.tolist()})'
761def _projection(vectors, direction):
762 dotprods = vectors @ direction
763 projection = direction[None, :] * dotprods[:, None]
764 return projection
767class FixedLine(IndexedConstraint):
768 """
769 Constrain an atom index or a list of atom indices to move on a line only.
771 The line is defined by its vector *direction*
772 """
774 def __init__(self, indices, direction):
775 """Constrain chosen atoms.
777 Parameters
778 ----------
779 indices : int or list of int
780 Index or indices for atoms that should be constrained
781 direction : list of 3 int
782 Direction of the vector defining the line
784 Examples
785 --------
786 Fix all Copper atoms to only move in the x-direction:
788 >>> from ase.constraints import FixedLine
789 >>> c = FixedLine(
790 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
791 ... direction=[1, 0, 0],
792 ... )
793 >>> atoms.set_constraint(c)
795 or constrain a single atom with the index 0 to move in the z-direction:
797 >>> c = FixedLine(indices=0, direction=[0, 0, 1])
798 >>> atoms.set_constraint(c)
799 """
800 super().__init__(indices)
801 self.dir = _normalize(direction)
803 def adjust_positions(self, atoms, newpositions):
804 step = newpositions[self.index] - atoms.positions[self.index]
805 projection = _projection(step, self.dir)
806 newpositions[self.index] = atoms.positions[self.index] + projection
808 def adjust_forces(self, atoms, forces):
809 forces[self.index] = _projection(forces[self.index], self.dir)
811 def get_removed_dof(self, atoms):
812 return 2 * len(self.index)
814 def __repr__(self):
815 return f'FixedLine(indices={self.index}, {self.dir.tolist()})'
817 def todict(self):
818 return {
819 'name': 'FixedLine',
820 'kwargs': {'indices': self.index.tolist(),
821 'direction': self.dir.tolist()}
822 }
825class FixCartesian(IndexedConstraint):
826 """Fix atoms in the directions of the cartesian coordinates.
828 Parameters
829 ----------
830 a : Sequence[int]
831 Indices of atoms to be fixed.
832 mask : tuple[bool, bool, bool], default: (True, True, True)
833 Cartesian directions to be fixed. (False: unfixed, True: fixed)
834 """
836 def __init__(self, a, mask=(True, True, True)):
837 super().__init__(indices=a)
838 self.mask = np.asarray(mask, bool)
840 def get_removed_dof(self, atoms: Atoms):
841 return self.mask.sum() * len(self.index)
843 def adjust_positions(self, atoms: Atoms, new):
844 new[self.index] = np.where(
845 self.mask[None, :],
846 atoms.positions[self.index],
847 new[self.index],
848 )
850 def adjust_forces(self, atoms: Atoms, forces):
851 forces[self.index] *= ~self.mask[None, :]
853 def todict(self):
854 return {'name': 'FixCartesian',
855 'kwargs': {'a': self.index.tolist(),
856 'mask': self.mask.tolist()}}
858 def __repr__(self):
859 name = type(self).__name__
860 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'
863class FixScaled(IndexedConstraint):
864 """Fix atoms in the directions of the unit vectors.
866 Parameters
867 ----------
868 a : Sequence[int]
869 Indices of atoms to be fixed.
870 mask : tuple[bool, bool, bool], default: (True, True, True)
871 Cell directions to be fixed. (False: unfixed, True: fixed)
872 """
874 def __init__(self, a, mask=(True, True, True), cell=None):
875 # XXX The unused cell keyword is there for compatibility
876 # with old trajectory files.
877 super().__init__(indices=a)
878 self.mask = np.asarray(mask, bool)
880 def get_removed_dof(self, atoms: Atoms):
881 return self.mask.sum() * len(self.index)
883 def adjust_positions(self, atoms: Atoms, new):
884 cell = atoms.cell
885 scaled_old = cell.scaled_positions(atoms.positions[self.index])
886 scaled_new = cell.scaled_positions(new[self.index])
887 scaled_new[:, self.mask] = scaled_old[:, self.mask]
888 new[self.index] = cell.cartesian_positions(scaled_new)
890 def adjust_forces(self, atoms: Atoms, forces):
891 # Forces are covariant to the coordinate transformation,
892 # use the inverse transformations
893 cell = atoms.cell
894 scaled_forces = cell.cartesian_positions(forces[self.index])
895 scaled_forces *= -(self.mask - 1)
896 forces[self.index] = cell.scaled_positions(scaled_forces)
898 def todict(self):
899 return {'name': 'FixScaled',
900 'kwargs': {'a': self.index.tolist(),
901 'mask': self.mask.tolist()}}
903 def __repr__(self):
904 name = type(self).__name__
905 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'
908# TODO: Better interface might be to use dictionaries in place of very
909# nested lists/tuples
910class FixInternals(FixConstraint):
911 """Constraint object for fixing multiple internal coordinates.
913 Allows fixing bonds, angles, dihedrals as well as linear combinations
914 of bonds (bondcombos).
916 Please provide angular units in degrees using `angles_deg` and
917 `dihedrals_deg`.
918 Fixing planar angles is not supported at the moment.
919 """
921 def __init__(self, bonds=None, angles=None, dihedrals=None,
922 angles_deg=None, dihedrals_deg=None,
923 bondcombos=None,
924 mic=False, epsilon=1.e-7):
925 """
926 A constrained internal coordinate is defined as a nested list:
927 '[value, [atom indices]]'. The constraint is initialized with a list of
928 constrained internal coordinates, i.e. '[[value, [atom indices]], ...]'.
929 If 'value' is None, the current value of the coordinate is constrained.
931 Parameters
932 ----------
933 bonds: nested python list, optional
934 List with targetvalue and atom indices defining the fixed bonds,
935 i.e. [[targetvalue, [index0, index1]], ...]
937 angles_deg: nested python list, optional
938 List with targetvalue and atom indices defining the fixedangles,
939 i.e. [[targetvalue, [index0, index1, index3]], ...]
941 dihedrals_deg: nested python list, optional
942 List with targetvalue and atom indices defining the fixed dihedrals,
943 i.e. [[targetvalue, [index0, index1, index3]], ...]
945 bondcombos: nested python list, optional
946 List with targetvalue, atom indices and linear coefficient defining
947 the fixed linear combination of bonds,
948 i.e. [[targetvalue, [[index0, index1, coefficient_for_bond],
949 [index1, index2, coefficient_for_bond]]], ...]
951 mic: bool, optional, default: False
952 Minimum image convention.
954 epsilon: float, optional, default: 1e-7
955 Convergence criterion.
956 """
957 warn_msg = 'Please specify {} in degrees using the {} argument.'
958 if angles:
959 warn(warn_msg.format('angles', 'angle_deg'), FutureWarning)
960 angles = np.asarray(angles)
961 angles[:, 0] = angles[:, 0] / np.pi * 180
962 angles = angles.tolist()
963 else:
964 angles = angles_deg
965 if dihedrals:
966 warn(warn_msg.format('dihedrals', 'dihedrals_deg'), FutureWarning)
967 dihedrals = np.asarray(dihedrals)
968 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180
969 dihedrals = dihedrals.tolist()
970 else:
971 dihedrals = dihedrals_deg
973 self.bonds = bonds or []
974 self.angles = angles or []
975 self.dihedrals = dihedrals or []
976 self.bondcombos = bondcombos or []
977 self.mic = mic
978 self.epsilon = epsilon
980 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)
981 + len(self.bondcombos))
983 # Initialize these at run-time:
984 self.constraints = []
985 self.initialized = False
987 def get_removed_dof(self, atoms):
988 return self.n
990 def initialize(self, atoms):
991 if self.initialized:
992 return
993 masses = np.repeat(atoms.get_masses(), 3)
994 cell = None
995 pbc = None
996 if self.mic:
997 cell = atoms.cell
998 pbc = atoms.pbc
999 self.constraints = []
1000 for data, ConstrClass in [(self.bonds, self.FixBondLengthAlt),
1001 (self.angles, self.FixAngle),
1002 (self.dihedrals, self.FixDihedral),
1003 (self.bondcombos, self.FixBondCombo)]:
1004 for datum in data:
1005 targetvalue = datum[0]
1006 if targetvalue is None: # set to current value
1007 targetvalue = ConstrClass.get_value(atoms, datum[1],
1008 self.mic)
1009 constr = ConstrClass(targetvalue, datum[1], masses, cell, pbc)
1010 self.constraints.append(constr)
1011 self.initialized = True
1013 @staticmethod
1014 def get_bondcombo(atoms, indices, mic=False):
1015 """Convenience function to return the value of the bondcombo coordinate
1016 (linear combination of bond lengths) for the given Atoms object 'atoms'.
1017 Example: Get the current value of the linear combination of two bond
1018 lengths defined as `bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]`."""
1019 c = sum(df[2] * atoms.get_distance(*df[:2], mic=mic) for df in indices)
1020 return c
1022 def get_subconstraint(self, atoms, definition):
1023 """Get pointer to a specific subconstraint.
1024 Identification by its definition via indices (and coefficients)."""
1025 self.initialize(atoms)
1026 for subconstr in self.constraints:
1027 if isinstance(definition[0], Sequence): # Combo constraint
1028 defin = [d + [c] for d, c in zip(subconstr.indices,
1029 subconstr.coefs)]
1030 if defin == definition:
1031 return subconstr
1032 else: # identify primitive constraints by their indices
1033 if subconstr.indices == [definition]:
1034 return subconstr
1035 raise ValueError('Given `definition` not found on Atoms object.')
1037 def shuffle_definitions(self, shuffle_dic, internal_type):
1038 dfns = [] # definitions
1039 for dfn in internal_type: # e.g. for bond in self.bonds
1040 append = True
1041 new_dfn = [dfn[0], list(dfn[1])]
1042 for old in dfn[1]:
1043 if old in shuffle_dic:
1044 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]
1045 else:
1046 append = False
1047 break
1048 if append:
1049 dfns.append(new_dfn)
1050 return dfns
1052 def shuffle_combos(self, shuffle_dic, internal_type):
1053 dfns = [] # definitions
1054 for dfn in internal_type: # i.e. for bondcombo in self.bondcombos
1055 append = True
1056 all_indices = [idx[0:-1] for idx in dfn[1]]
1057 new_dfn = [dfn[0], list(dfn[1])]
1058 for i, indices in enumerate(all_indices):
1059 for old in indices:
1060 if old in shuffle_dic:
1061 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]
1062 else:
1063 append = False
1064 break
1065 if not append:
1066 break
1067 if append:
1068 dfns.append(new_dfn)
1069 return dfns
1071 def index_shuffle(self, atoms, ind):
1072 # See docstring of superclass
1073 self.initialize(atoms)
1074 shuffle_dic = dict(slice2enlist(ind, len(atoms)))
1075 shuffle_dic = {old: new for new, old in shuffle_dic.items()}
1076 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)
1077 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)
1078 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)
1079 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)
1080 self.initialized = False
1081 self.initialize(atoms)
1082 if len(self.constraints) == 0:
1083 raise IndexError('Constraint not part of slice')
1085 def get_indices(self):
1086 cons = []
1087 for dfn in self.bonds + self.dihedrals + self.angles:
1088 cons.extend(dfn[1])
1089 for dfn in self.bondcombos:
1090 for partial_dfn in dfn[1]:
1091 cons.extend(partial_dfn[0:-1]) # last index is the coefficient
1092 return list(set(cons))
1094 def todict(self):
1095 return {'name': 'FixInternals',
1096 'kwargs': {'bonds': self.bonds,
1097 'angles_deg': self.angles,
1098 'dihedrals_deg': self.dihedrals,
1099 'bondcombos': self.bondcombos,
1100 'mic': self.mic,
1101 'epsilon': self.epsilon}}
1103 def adjust_positions(self, atoms, newpos):
1104 self.initialize(atoms)
1105 for constraint in self.constraints:
1106 constraint.setup_jacobian(atoms.positions)
1107 for _ in range(50):
1108 maxerr = 0.0
1109 for constraint in self.constraints:
1110 constraint.adjust_positions(atoms.positions, newpos)
1111 maxerr = max(abs(constraint.sigma), maxerr)
1112 if maxerr < self.epsilon:
1113 return
1114 msg = 'FixInternals.adjust_positions did not converge.'
1115 if any(constr.targetvalue > 175. or constr.targetvalue < 5. for constr
1116 in self.constraints if isinstance(constr, self.FixAngle)):
1117 msg += (' This may be caused by an almost planar angle.'
1118 ' Support for planar angles would require the'
1119 ' implementation of ghost, i.e. dummy, atoms.'
1120 ' See issue #868.')
1121 raise ValueError(msg)
1123 def adjust_forces(self, atoms, forces):
1124 """Project out translations and rotations and all other constraints"""
1125 self.initialize(atoms)
1126 positions = atoms.positions
1127 N = len(forces)
1128 list2_constraints = list(np.zeros((6, N, 3)))
1129 tx, ty, tz, rx, ry, rz = list2_constraints
1131 list_constraints = [r.ravel() for r in list2_constraints]
1133 tx[:, 0] = 1.0
1134 ty[:, 1] = 1.0
1135 tz[:, 2] = 1.0
1136 ff = forces.ravel()
1138 # Calculate the center of mass
1139 center = positions.sum(axis=0) / N
1141 rx[:, 1] = -(positions[:, 2] - center[2])
1142 rx[:, 2] = positions[:, 1] - center[1]
1143 ry[:, 0] = positions[:, 2] - center[2]
1144 ry[:, 2] = -(positions[:, 0] - center[0])
1145 rz[:, 0] = -(positions[:, 1] - center[1])
1146 rz[:, 1] = positions[:, 0] - center[0]
1148 # Normalizing transl., rotat. constraints
1149 for r in list2_constraints:
1150 r /= np.linalg.norm(r.ravel())
1152 # Add all angle, etc. constraint vectors
1153 for constraint in self.constraints:
1154 constraint.setup_jacobian(positions)
1155 constraint.adjust_forces(positions, forces)
1156 list_constraints.insert(0, constraint.jacobian)
1157 # QR DECOMPOSITION - GRAM SCHMIDT
1159 list_constraints = [r.ravel() for r in list_constraints]
1160 aa = np.column_stack(list_constraints)
1161 (aa, bb) = np.linalg.qr(aa)
1162 # Projection
1163 hh = []
1164 for i, constraint in enumerate(self.constraints):
1165 hh.append(aa[:, i] * np.vstack(aa[:, i]))
1167 txx = aa[:, self.n] * np.vstack(aa[:, self.n])
1168 tyy = aa[:, self.n + 1] * np.vstack(aa[:, self.n + 1])
1169 tzz = aa[:, self.n + 2] * np.vstack(aa[:, self.n + 2])
1170 rxx = aa[:, self.n + 3] * np.vstack(aa[:, self.n + 3])
1171 ryy = aa[:, self.n + 4] * np.vstack(aa[:, self.n + 4])
1172 rzz = aa[:, self.n + 5] * np.vstack(aa[:, self.n + 5])
1173 T = txx + tyy + tzz + rxx + ryy + rzz
1174 for vec in hh:
1175 T += vec
1176 ff = np.dot(T, np.vstack(ff))
1177 forces[:, :] -= np.dot(T, np.vstack(ff)).reshape(-1, 3)
1179 def __repr__(self):
1180 constraints = [repr(constr) for constr in self.constraints]
1181 return f'FixInternals(_copy_init={constraints}, epsilon={self.epsilon})'
1183 # Classes for internal use in FixInternals
1184 class FixInternalsBase:
1185 """Base class for subclasses of FixInternals."""
1187 def __init__(self, targetvalue, indices, masses, cell, pbc):
1188 self.targetvalue = targetvalue # constant target value
1189 self.indices = [defin[0:-1] for defin in indices] # indices, defs
1190 self.coefs = np.asarray([defin[-1] for defin in indices])
1191 self.masses = masses
1192 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix
1193 self.sigma = 1. # difference between current and target value
1194 self.projected_force = None # helps optimizers scan along constr.
1195 self.cell = cell
1196 self.pbc = pbc
1198 def finalize_jacobian(self, pos, n_internals, n, derivs):
1199 """Populate jacobian with derivatives for `n_internals` defined
1200 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""
1201 jacobian = np.zeros((n_internals, *pos.shape))
1202 for i, idx in enumerate(self.indices):
1203 for j in range(n):
1204 jacobian[i, idx[j]] = derivs[i, j]
1205 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))
1206 return self.coefs @ jacobian
1208 def finalize_positions(self, newpos):
1209 jacobian = self.jacobian / self.masses
1210 lamda = -self.sigma / (jacobian @ self.get_jacobian(newpos))
1211 dnewpos = lamda * jacobian
1212 newpos += dnewpos.reshape(newpos.shape)
1214 def adjust_forces(self, positions, forces):
1215 self.projected_forces = ((self.jacobian @ forces.ravel())
1216 * self.jacobian)
1217 self.jacobian /= np.linalg.norm(self.jacobian)
1219 class FixBondCombo(FixInternalsBase):
1220 """Constraint subobject for fixing linear combination of bond lengths
1221 within FixInternals.
1223 sum_i( coef_i * bond_length_i ) = constant
1224 """
1226 def get_jacobian(self, pos):
1227 bondvectors = [pos[k] - pos[h] for h, k in self.indices]
1228 derivs = get_distances_derivatives(bondvectors, cell=self.cell,
1229 pbc=self.pbc)
1230 return self.finalize_jacobian(pos, len(bondvectors), 2, derivs)
1232 def setup_jacobian(self, pos):
1233 self.jacobian = self.get_jacobian(pos)
1235 def adjust_positions(self, oldpos, newpos):
1236 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]
1237 (_, ), (dists, ) = conditional_find_mic([bondvectors],
1238 cell=self.cell,
1239 pbc=self.pbc)
1240 value = self.coefs @ dists
1241 self.sigma = value - self.targetvalue
1242 self.finalize_positions(newpos)
1244 @staticmethod
1245 def get_value(atoms, indices, mic):
1246 return FixInternals.get_bondcombo(atoms, indices, mic)
1248 def __repr__(self):
1249 return (f'FixBondCombo({self.targetvalue}, {self.indices}, '
1250 '{self.coefs})')
1252 class FixBondLengthAlt(FixBondCombo):
1253 """Constraint subobject for fixing bond length within FixInternals.
1254 Fix distance between atoms with indices a1, a2."""
1256 def __init__(self, targetvalue, indices, masses, cell, pbc):
1257 if targetvalue <= 0.:
1258 raise ZeroDivisionError('Invalid targetvalue for fixed bond')
1259 indices = [list(indices) + [1.]] # bond definition with coef 1.
1260 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1262 @staticmethod
1263 def get_value(atoms, indices, mic):
1264 return atoms.get_distance(*indices, mic=mic)
1266 def __repr__(self):
1267 return f'FixBondLengthAlt({self.targetvalue}, {self.indices})'
1269 class FixAngle(FixInternalsBase):
1270 """Constraint subobject for fixing an angle within FixInternals.
1272 Convergence is potentially problematic for angles very close to
1273 0 or 180 degrees as there is a singularity in the Cartesian derivative.
1274 Fixing planar angles is therefore not supported at the moment.
1275 """
1277 def __init__(self, targetvalue, indices, masses, cell, pbc):
1278 """Fix atom movement to construct a constant angle."""
1279 if targetvalue <= 0. or targetvalue >= 180.:
1280 raise ZeroDivisionError('Invalid targetvalue for fixed angle')
1281 indices = [list(indices) + [1.]] # angle definition with coef 1.
1282 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1284 def gather_vectors(self, pos):
1285 v0 = [pos[h] - pos[k] for h, k, l in self.indices]
1286 v1 = [pos[l] - pos[k] for h, k, l in self.indices]
1287 return v0, v1
1289 def get_jacobian(self, pos):
1290 v0, v1 = self.gather_vectors(pos)
1291 derivs = get_angles_derivatives(v0, v1, cell=self.cell,
1292 pbc=self.pbc)
1293 return self.finalize_jacobian(pos, len(v0), 3, derivs)
1295 def setup_jacobian(self, pos):
1296 self.jacobian = self.get_jacobian(pos)
1298 def adjust_positions(self, oldpos, newpos):
1299 v0, v1 = self.gather_vectors(newpos)
1300 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)
1301 self.sigma = value - self.targetvalue
1302 self.finalize_positions(newpos)
1304 @staticmethod
1305 def get_value(atoms, indices, mic):
1306 return atoms.get_angle(*indices, mic=mic)
1308 def __repr__(self):
1309 return f'FixAngle({self.targetvalue}, {self.indices})'
1311 class FixDihedral(FixInternalsBase):
1312 """Constraint subobject for fixing a dihedral angle within FixInternals.
1314 A dihedral becomes undefined when at least one of the inner two angles
1315 becomes planar. Make sure to avoid this situation.
1316 """
1318 def __init__(self, targetvalue, indices, masses, cell, pbc):
1319 indices = [list(indices) + [1.]] # dihedral def. with coef 1.
1320 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1322 def gather_vectors(self, pos):
1323 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]
1324 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]
1325 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]
1326 return v0, v1, v2
1328 def get_jacobian(self, pos):
1329 v0, v1, v2 = self.gather_vectors(pos)
1330 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,
1331 pbc=self.pbc)
1332 return self.finalize_jacobian(pos, len(v0), 4, derivs)
1334 def setup_jacobian(self, pos):
1335 self.jacobian = self.get_jacobian(pos)
1337 def adjust_positions(self, oldpos, newpos):
1338 v0, v1, v2 = self.gather_vectors(newpos)
1339 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1340 # apply minimum dihedral difference 'convention': (diff <= 180)
1341 self.sigma = (value - self.targetvalue + 180) % 360 - 180
1342 self.finalize_positions(newpos)
1344 @staticmethod
1345 def get_value(atoms, indices, mic):
1346 return atoms.get_dihedral(*indices, mic=mic)
1348 def __repr__(self):
1349 return f'FixDihedral({self.targetvalue}, {self.indices})'
1352class FixParametricRelations(FixConstraint):
1354 def __init__(
1355 self,
1356 indices,
1357 Jacobian,
1358 const_shift,
1359 params=None,
1360 eps=1e-12,
1361 use_cell=False,
1362 ):
1363 """Constrains the degrees of freedom to act in a reduced parameter
1364 space defined by the Jacobian
1366 These constraints are based off the work in:
1367 https://arxiv.org/abs/1908.01610
1369 The constraints linearly maps the full 3N degrees of freedom,
1370 where N is number of active lattice vectors/atoms onto a
1371 reduced subset of M free parameters, where M <= 3*N. The
1372 Jacobian matrix and constant shift vector map the full set of
1373 degrees of freedom onto the reduced parameter space.
1375 Currently the constraint is set up to handle either atomic
1376 positions or lattice vectors at one time, but not both. To do
1377 both simply add a two constraints for each set. This is done
1378 to keep the mathematics behind the operations separate.
1380 It would be possible to extend these constraints to allow
1381 non-linear transformations if functionality to update the
1382 Jacobian at each position update was included. This would
1383 require passing an update function evaluate it every time
1384 adjust_positions is callled. This is currently NOT supported,
1385 and there are no plans to implement it in the future.
1387 Args:
1388 indices (list of int): indices of the constrained atoms
1389 (if not None or empty then cell_indices must be None or Empty)
1390 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))):
1391 The Jacobian describing
1392 the parameter space transformation
1393 const_shift (np.ndarray(shape=(3*len(indices)))):
1394 A vector describing the constant term
1395 in the transformation not accounted for in the Jacobian
1396 params (list of str):
1397 parameters used in the parametric representation
1398 if None a list is generated based on the shape of the Jacobian
1399 eps (float): a small number to compare the similarity of
1400 numbers and set the precision used
1401 to generate the constraint expressions
1402 use_cell (bool): if True then act on the cell object
1404 """
1405 self.indices = np.array(indices)
1406 self.Jacobian = np.array(Jacobian)
1407 self.const_shift = np.array(const_shift)
1409 assert self.const_shift.shape[0] == 3 * len(self.indices)
1410 assert self.Jacobian.shape[0] == 3 * len(self.indices)
1412 self.eps = eps
1413 self.use_cell = use_cell
1415 if params is None:
1416 params = []
1417 if self.Jacobian.shape[1] > 0:
1418 int_fmt_str = "{:0" + \
1419 str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"
1420 for param_ind in range(self.Jacobian.shape[1]):
1421 params.append("param_" + int_fmt_str.format(param_ind))
1422 else:
1423 assert len(params) == self.Jacobian.shape[-1]
1425 self.params = params
1427 self.Jacobian_inv = np.linalg.inv(
1428 self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T
1430 @classmethod
1431 def from_expressions(cls, indices, params, expressions,
1432 eps=1e-12, use_cell=False):
1433 """Converts the expressions into a Jacobian Matrix/const_shift
1434 vector and constructs a FixParametricRelations constraint
1436 The expressions must be a list like object of size 3*N and
1437 elements must be ordered as:
1438 [n_0,i; n_0,j; n_0,k; n_1,i; n_1,j; .... ; n_N-1,i; n_N-1,j; n_N-1,k],
1439 where i, j, and k are the first, second and third
1440 component of the atomic position/lattice
1441 vector. Currently only linear operations are allowed to be
1442 included in the expressions so
1443 only terms like:
1444 - const * param_0
1445 - sqrt[const] * param_1
1446 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M
1447 where const is any real number and param_0, param_1, ..., param_M are
1448 the parameters passed in
1449 params, are allowed.
1451 For example, fractional atomic position constraints for wurtzite are:
1452 params = ["z1", "z2"]
1453 expressions = [
1454 "1.0/3.0", "2.0/3.0", "z1",
1455 "2.0/3.0", "1.0/3.0", "0.5 + z1",
1456 "1.0/3.0", "2.0/3.0", "z2",
1457 "2.0/3.0", "1.0/3.0", "0.5 + z2",
1458 ]
1460 For diamond are:
1461 params = []
1462 expressions = [
1463 "0.0", "0.0", "0.0",
1464 "0.25", "0.25", "0.25",
1465 ],
1467 and for stannite are
1468 params=["x4", "z4"]
1469 expressions = [
1470 "0.0", "0.0", "0.0",
1471 "0.0", "0.5", "0.5",
1472 "0.75", "0.25", "0.5",
1473 "0.25", "0.75", "0.5",
1474 "x4 + z4", "x4 + z4", "2*x4",
1475 "x4 - z4", "x4 - z4", "-2*x4",
1476 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",
1477 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",
1478 ]
1480 Args:
1481 indices (list of int): indices of the constrained atoms
1482 (if not None or empty then cell_indices must be None or Empty)
1483 params (list of str): parameters used in the
1484 parametric representation
1485 expressions (list of str): expressions used to convert from the
1486 parametric to the real space representation
1487 eps (float): a small number to compare the similarity of
1488 numbers and set the precision used
1489 to generate the constraint expressions
1490 use_cell (bool): if True then act on the cell object
1492 Returns:
1493 cls(
1494 indices,
1495 Jacobian generated from expressions,
1496 const_shift generated from expressions,
1497 params,
1498 eps-12,
1499 use_cell,
1500 )
1501 """
1502 Jacobian = np.zeros((3 * len(indices), len(params)))
1503 const_shift = np.zeros(3 * len(indices))
1505 for expr_ind, expression in enumerate(expressions):
1506 expression = expression.strip()
1508 # Convert subtraction to addition
1509 expression = expression.replace("-", "+(-1.0)*")
1510 if expression[0] == "+":
1511 expression = expression[1:]
1512 elif expression[:2] == "(+":
1513 expression = "(" + expression[2:]
1515 # Explicitly add leading zeros so when replacing param_1 with 0.0
1516 # param_11 does not become 0.01
1517 int_fmt_str = "{:0" + \
1518 str(int(np.ceil(np.log10(len(params) + 1)))) + "d}"
1520 param_dct = {}
1521 param_map = {}
1523 # Construct a standardized param template for A/B filling
1524 for param_ind, param in enumerate(params):
1525 param_str = "param_" + int_fmt_str.format(param_ind)
1526 param_map[param] = param_str
1527 param_dct[param_str] = 0.0
1529 # Replace the parameters according to the map
1530 # Sort by string length (long to short) to prevent cases like x11
1531 # becoming f"{param_map["x1"]}1"
1532 for param in sorted(params, key=lambda s: -1.0 * len(s)):
1533 expression = expression.replace(param, param_map[param])
1535 # Partial linearity check
1536 for express_sec in expression.split("+"):
1537 in_sec = [param in express_sec for param in param_dct]
1538 n_params_in_sec = len(np.where(np.array(in_sec))[0])
1539 if n_params_in_sec > 1:
1540 raise ValueError(
1541 "FixParametricRelations expressions must be linear.")
1543 const_shift[expr_ind] = float(
1544 eval_expression(expression, param_dct))
1546 for param_ind in range(len(params)):
1547 param_str = "param_" + int_fmt_str.format(param_ind)
1548 if param_str not in expression:
1549 Jacobian[expr_ind, param_ind] = 0.0
1550 continue
1551 param_dct[param_str] = 1.0
1552 test_1 = float(eval_expression(expression, param_dct))
1553 test_1 -= const_shift[expr_ind]
1554 Jacobian[expr_ind, param_ind] = test_1
1556 param_dct[param_str] = 2.0
1557 test_2 = float(eval_expression(expression, param_dct))
1558 test_2 -= const_shift[expr_ind]
1559 if abs(test_2 / test_1 - 2.0) > eps:
1560 raise ValueError(
1561 "FixParametricRelations expressions must be linear.")
1562 param_dct[param_str] = 0.0
1564 args = [
1565 indices,
1566 Jacobian,
1567 const_shift,
1568 params,
1569 eps,
1570 use_cell,
1571 ]
1572 if cls is FixScaledParametricRelations:
1573 args = args[:-1]
1574 return cls(*args)
1576 @property
1577 def expressions(self):
1578 """Generate the expressions represented by the current self.Jacobian
1579 and self.const_shift objects"""
1580 expressions = []
1581 per = int(round(-1 * np.log10(self.eps)))
1582 fmt_str = "{:." + str(per + 1) + "g}"
1583 for index, shift_val in enumerate(self.const_shift):
1584 exp = ""
1585 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(
1586 shift_val) > self.eps:
1587 exp += fmt_str.format(shift_val)
1589 param_exp = ""
1590 for param_index, jacob_val in enumerate(self.Jacobian[index]):
1591 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)
1592 if abs_jacob_val < self.eps:
1593 continue
1595 param = self.params[param_index]
1596 if param_exp or exp:
1597 if jacob_val > -1.0 * self.eps:
1598 param_exp += " + "
1599 else:
1600 param_exp += " - "
1601 elif (not exp) and (not param_exp) and (
1602 jacob_val < -1.0 * self.eps):
1603 param_exp += "-"
1605 if np.abs(abs_jacob_val - 1.0) <= self.eps:
1606 param_exp += f"{param:s}"
1607 else:
1608 param_exp += (fmt_str +
1609 "*{:s}").format(abs_jacob_val, param)
1611 exp += param_exp
1613 expressions.append(exp)
1614 return np.array(expressions).reshape((-1, 3))
1616 def todict(self):
1617 """Create a dictionary representation of the constraint"""
1618 return {
1619 "name": type(self).__name__,
1620 "kwargs": {
1621 "indices": self.indices,
1622 "params": self.params,
1623 "Jacobian": self.Jacobian,
1624 "const_shift": self.const_shift,
1625 "eps": self.eps,
1626 "use_cell": self.use_cell,
1627 }
1628 }
1630 def __repr__(self):
1631 """The str representation of the constraint"""
1632 if len(self.indices) > 1:
1633 indices_str = "[{:d}, ..., {:d}]".format(
1634 self.indices[0], self.indices[-1])
1635 else:
1636 indices_str = f"[{self.indices[0]:d}]"
1638 if len(self.params) > 1:
1639 params_str = "[{:s}, ..., {:s}]".format(
1640 self.params[0], self.params[-1])
1641 elif len(self.params) == 1:
1642 params_str = f"[{self.params[0]:s}]"
1643 else:
1644 params_str = "[]"
1646 return '{:s}({:s}, {:s}, ..., {:e})'.format(
1647 type(self).__name__,
1648 indices_str,
1649 params_str,
1650 self.eps
1651 )
1654class FixScaledParametricRelations(FixParametricRelations):
1656 def __init__(
1657 self,
1658 indices,
1659 Jacobian,
1660 const_shift,
1661 params=None,
1662 eps=1e-12,
1663 ):
1664 """The fractional coordinate version of FixParametricRelations
1666 All arguments are the same, but since this is for fractional
1667 coordinates use_cell is false"""
1668 super().__init__(
1669 indices,
1670 Jacobian,
1671 const_shift,
1672 params,
1673 eps,
1674 False,
1675 )
1677 def adjust_contravariant(self, cell, vecs, B):
1678 """Adjust the values of a set of vectors that are contravariant
1679 with the unit transformation"""
1680 scaled = cell.scaled_positions(vecs).flatten()
1681 scaled = self.Jacobian_inv @ (scaled - B)
1682 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))
1684 return cell.cartesian_positions(scaled)
1686 def adjust_positions(self, atoms, positions):
1687 """Adjust positions of the atoms to match the constraints"""
1688 positions[self.indices] = self.adjust_contravariant(
1689 atoms.cell,
1690 positions[self.indices],
1691 self.const_shift,
1692 )
1693 positions[self.indices] = self.adjust_B(
1694 atoms.cell, positions[self.indices])
1696 def adjust_B(self, cell, positions):
1697 """Wraps the positions back to the unit cell and adjust B to
1698 keep track of this change"""
1699 fractional = cell.scaled_positions(positions)
1700 wrapped_fractional = (fractional % 1.0) % 1.0
1701 self.const_shift += np.round(wrapped_fractional - fractional).flatten()
1702 return cell.cartesian_positions(wrapped_fractional)
1704 def adjust_momenta(self, atoms, momenta):
1705 """Adjust momenta of the atoms to match the constraints"""
1706 momenta[self.indices] = self.adjust_contravariant(
1707 atoms.cell,
1708 momenta[self.indices],
1709 np.zeros(self.const_shift.shape),
1710 )
1712 def adjust_forces(self, atoms, forces):
1713 """Adjust forces of the atoms to match the constraints"""
1714 # Forces are coavarient to the coordinate transformation, use the
1715 # inverse transformations
1716 cart2frac_jacob = np.zeros(2 * (3 * len(atoms),))
1717 for i_atom in range(len(atoms)):
1718 cart2frac_jacob[3 * i_atom:3 * (i_atom + 1),
1719 3 * i_atom:3 * (i_atom + 1)] = atoms.cell.T
1721 jacobian = cart2frac_jacob @ self.Jacobian
1722 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T
1724 reduced_forces = jacobian.T @ forces.flatten()
1725 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)
1727 def todict(self):
1728 """Create a dictionary representation of the constraint"""
1729 dct = super().todict()
1730 del dct["kwargs"]["use_cell"]
1731 return dct
1734class FixCartesianParametricRelations(FixParametricRelations):
1736 def __init__(
1737 self,
1738 indices,
1739 Jacobian,
1740 const_shift,
1741 params=None,
1742 eps=1e-12,
1743 use_cell=False,
1744 ):
1745 """The Cartesian coordinate version of FixParametricRelations"""
1746 super().__init__(
1747 indices,
1748 Jacobian,
1749 const_shift,
1750 params,
1751 eps,
1752 use_cell,
1753 )
1755 def adjust_contravariant(self, vecs, B):
1756 """Adjust the values of a set of vectors that are contravariant with
1757 the unit transformation"""
1758 vecs = self.Jacobian_inv @ (vecs.flatten() - B)
1759 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))
1760 return vecs
1762 def adjust_positions(self, atoms, positions):
1763 """Adjust positions of the atoms to match the constraints"""
1764 if self.use_cell:
1765 return
1766 positions[self.indices] = self.adjust_contravariant(
1767 positions[self.indices],
1768 self.const_shift,
1769 )
1771 def adjust_momenta(self, atoms, momenta):
1772 """Adjust momenta of the atoms to match the constraints"""
1773 if self.use_cell:
1774 return
1775 momenta[self.indices] = self.adjust_contravariant(
1776 momenta[self.indices],
1777 np.zeros(self.const_shift.shape),
1778 )
1780 def adjust_forces(self, atoms, forces):
1781 """Adjust forces of the atoms to match the constraints"""
1782 if self.use_cell:
1783 return
1785 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()
1786 forces[self.indices] = (self.Jacobian_inv.T @
1787 forces_reduced).reshape(-1, 3)
1789 def adjust_cell(self, atoms, cell):
1790 """Adjust the cell of the atoms to match the constraints"""
1791 if not self.use_cell:
1792 return
1793 cell[self.indices] = self.adjust_contravariant(
1794 cell[self.indices],
1795 np.zeros(self.const_shift.shape),
1796 )
1798 def adjust_stress(self, atoms, stress):
1799 """Adjust the stress of the atoms to match the constraints"""
1800 if not self.use_cell:
1801 return
1803 stress_3x3 = voigt_6_to_full_3x3_stress(stress)
1804 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()
1805 stress_3x3[self.indices] = (
1806 self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)
1808 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)
1811class Hookean(FixConstraint):
1812 """Applies a Hookean restorative force between a pair of atoms, an atom
1813 and a point, or an atom and a plane."""
1815 def __init__(self, a1, a2, k, rt=None):
1816 """Forces two atoms to stay close together by applying no force if
1817 they are below a threshold length, rt, and applying a Hookean
1818 restorative force when the distance between them exceeds rt. Can
1819 also be used to tether an atom to a fixed point in space or to a
1820 distance above a plane.
1822 a1 : int
1823 Index of atom 1
1824 a2 : one of three options
1825 1) index of atom 2
1826 2) a fixed point in cartesian space to which to tether a1
1827 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.
1828 k : float
1829 Hooke's law (spring) constant to apply when distance
1830 exceeds threshold_length. Units of eV A^-2.
1831 rt : float
1832 The threshold length below which there is no force. The
1833 length is 1) between two atoms, 2) between atom and point.
1834 This argument is not supplied in case 3. Units of A.
1836 If a plane is specified, the Hooke's law force is applied if the atom
1837 is on the normal side of the plane. For instance, the plane with
1838 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z
1839 intercept of +7 and a normal vector pointing in the +z direction.
1840 If the atom has z > 7, then a downward force would be applied of
1841 k * (atom.z - 7). The same plane with the normal vector pointing in
1842 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).
1844 References:
1846 Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014)
1847 https://link.springer.com/article/10.1007%2Fs11244-013-0161-8
1848 """
1850 if isinstance(a2, int):
1851 self._type = 'two atoms'
1852 self.indices = [a1, a2]
1853 elif len(a2) == 3:
1854 self._type = 'point'
1855 self.index = a1
1856 self.origin = np.array(a2)
1857 elif len(a2) == 4:
1858 self._type = 'plane'
1859 self.index = a1
1860 self.plane = a2
1861 else:
1862 raise RuntimeError('Unknown type for a2')
1863 self.threshold = rt
1864 self.spring = k
1866 def get_removed_dof(self, atoms):
1867 return 0
1869 def todict(self):
1870 dct = {'name': 'Hookean'}
1871 dct['kwargs'] = {'rt': self.threshold,
1872 'k': self.spring}
1873 if self._type == 'two atoms':
1874 dct['kwargs']['a1'] = self.indices[0]
1875 dct['kwargs']['a2'] = self.indices[1]
1876 elif self._type == 'point':
1877 dct['kwargs']['a1'] = self.index
1878 dct['kwargs']['a2'] = self.origin
1879 elif self._type == 'plane':
1880 dct['kwargs']['a1'] = self.index
1881 dct['kwargs']['a2'] = self.plane
1882 else:
1883 raise NotImplementedError(f'Bad type: {self._type}')
1884 return dct
1886 def adjust_positions(self, atoms, newpositions):
1887 pass
1889 def adjust_momenta(self, atoms, momenta):
1890 pass
1892 def adjust_forces(self, atoms, forces):
1893 positions = atoms.positions
1894 if self._type == 'plane':
1895 A, B, C, D = self.plane
1896 x, y, z = positions[self.index]
1897 d = ((A * x + B * y + C * z + D) /
1898 np.sqrt(A**2 + B**2 + C**2))
1899 if d < 0:
1900 return
1901 magnitude = self.spring * d
1902 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))
1903 forces[self.index] += direction * magnitude
1904 return
1905 if self._type == 'two atoms':
1906 p1, p2 = positions[self.indices]
1907 elif self._type == 'point':
1908 p1 = positions[self.index]
1909 p2 = self.origin
1910 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1911 bondlength = np.linalg.norm(displace)
1912 if bondlength > self.threshold:
1913 magnitude = self.spring * (bondlength - self.threshold)
1914 direction = displace / np.linalg.norm(displace)
1915 if self._type == 'two atoms':
1916 forces[self.indices[0]] += direction * magnitude
1917 forces[self.indices[1]] -= direction * magnitude
1918 else:
1919 forces[self.index] += direction * magnitude
1921 def adjust_potential_energy(self, atoms):
1922 """Returns the difference to the potential energy due to an active
1923 constraint. (That is, the quantity returned is to be added to the
1924 potential energy.)"""
1925 positions = atoms.positions
1926 if self._type == 'plane':
1927 A, B, C, D = self.plane
1928 x, y, z = positions[self.index]
1929 d = ((A * x + B * y + C * z + D) /
1930 np.sqrt(A**2 + B**2 + C**2))
1931 if d > 0:
1932 return 0.5 * self.spring * d**2
1933 else:
1934 return 0.
1935 if self._type == 'two atoms':
1936 p1, p2 = positions[self.indices]
1937 elif self._type == 'point':
1938 p1 = positions[self.index]
1939 p2 = self.origin
1940 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1941 bondlength = np.linalg.norm(displace)
1942 if bondlength > self.threshold:
1943 return 0.5 * self.spring * (bondlength - self.threshold)**2
1944 else:
1945 return 0.
1947 def get_indices(self):
1948 if self._type == 'two atoms':
1949 return self.indices
1950 elif self._type == 'point':
1951 return self.index
1952 elif self._type == 'plane':
1953 return self.index
1955 def index_shuffle(self, atoms, ind):
1956 # See docstring of superclass
1957 if self._type == 'two atoms':
1958 newa = [-1, -1] # Signal error
1959 for new, old in slice2enlist(ind, len(atoms)):
1960 for i, a in enumerate(self.indices):
1961 if old == a:
1962 newa[i] = new
1963 if newa[0] == -1 or newa[1] == -1:
1964 raise IndexError('Constraint not part of slice')
1965 self.indices = newa
1966 elif (self._type == 'point') or (self._type == 'plane'):
1967 newa = -1 # Signal error
1968 for new, old in slice2enlist(ind, len(atoms)):
1969 if old == self.index:
1970 newa = new
1971 break
1972 if newa == -1:
1973 raise IndexError('Constraint not part of slice')
1974 self.index = newa
1976 def __repr__(self):
1977 if self._type == 'two atoms':
1978 return 'Hookean(%d, %d)' % tuple(self.indices)
1979 elif self._type == 'point':
1980 return 'Hookean(%d) to cartesian' % self.index
1981 else:
1982 return 'Hookean(%d) to plane' % self.index
1985class ExternalForce(FixConstraint):
1986 """Constraint object for pulling two atoms apart by an external force.
1988 You can combine this constraint for example with FixBondLength but make
1989 sure that *ExternalForce* comes first in the list if there are overlaps
1990 between atom1-2 and atom3-4:
1992 >>> from ase.build import bulk
1994 >>> atoms = bulk('Cu', 'fcc', a=3.6)
1995 >>> atom1, atom2, atom3, atom4 = atoms[:4]
1996 >>> fext = 1.0
1997 >>> con1 = ExternalForce(atom1, atom2, f_ext)
1998 >>> con2 = FixBondLength(atom3, atom4)
1999 >>> atoms.set_constraint([con1, con2])
2001 see ase/test/external_force.py"""
2003 def __init__(self, a1, a2, f_ext):
2004 self.indices = [a1, a2]
2005 self.external_force = f_ext
2007 def get_removed_dof(self, atoms):
2008 return 0
2010 def adjust_positions(self, atoms, new):
2011 pass
2013 def adjust_forces(self, atoms, forces):
2014 dist = np.subtract.reduce(atoms.positions[self.indices])
2015 force = self.external_force * dist / np.linalg.norm(dist)
2016 forces[self.indices] += (force, -force)
2018 def adjust_potential_energy(self, atoms):
2019 dist = np.subtract.reduce(atoms.positions[self.indices])
2020 return -np.linalg.norm(dist) * self.external_force
2022 def index_shuffle(self, atoms, ind):
2023 """Shuffle the indices of the two atoms in this constraint"""
2024 newa = [-1, -1] # Signal error
2025 for new, old in slice2enlist(ind, len(atoms)):
2026 for i, a in enumerate(self.indices):
2027 if old == a:
2028 newa[i] = new
2029 if newa[0] == -1 or newa[1] == -1:
2030 raise IndexError('Constraint not part of slice')
2031 self.indices = newa
2033 def __repr__(self):
2034 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],
2035 self.indices[1],
2036 self.external_force)
2038 def todict(self):
2039 return {'name': 'ExternalForce',
2040 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2041 'f_ext': self.external_force}}
2044class MirrorForce(FixConstraint):
2045 """Constraint object for mirroring the force between two atoms.
2047 This class is designed to find a transition state with the help of a
2048 single optimization. It can be used if the transition state belongs to a
2049 bond breaking reaction. First the given bond length will be fixed until
2050 all other degrees of freedom are optimized, then the forces of the two
2051 atoms will be mirrored to find the transition state. The mirror plane is
2052 perpendicular to the connecting line of the atoms. Transition states in
2053 dependence of the force can be obtained by stretching the molecule and
2054 fixing its total length with *FixBondLength* or by using *ExternalForce*
2055 during the optimization with *MirrorForce*.
2057 Parameters
2058 ----------
2059 a1: int
2060 First atom index.
2061 a2: int
2062 Second atom index.
2063 max_dist: float
2064 Upper limit of the bond length interval where the transition state
2065 can be found.
2066 min_dist: float
2067 Lower limit of the bond length interval where the transition state
2068 can be found.
2069 fmax: float
2070 Maximum force used for the optimization.
2072 Notes
2073 -----
2074 You can combine this constraint for example with FixBondLength but make
2075 sure that *MirrorForce* comes first in the list if there are overlaps
2076 between atom1-2 and atom3-4:
2078 >>> from ase.build import bulk
2080 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2081 >>> atom1, atom2, atom3, atom4 = atoms[:4]
2082 >>> con1 = MirrorForce(atom1, atom2)
2083 >>> con2 = FixBondLength(atom3, atom4)
2084 >>> atoms.set_constraint([con1, con2])
2086 """
2088 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):
2089 self.indices = [a1, a2]
2090 self.min_dist = min_dist
2091 self.max_dist = max_dist
2092 self.fmax = fmax
2094 def adjust_positions(self, atoms, new):
2095 pass
2097 def adjust_forces(self, atoms, forces):
2098 dist = np.subtract.reduce(atoms.positions[self.indices])
2099 d = np.linalg.norm(dist)
2100 if (d < self.min_dist) or (d > self.max_dist):
2101 # Stop structure optimization
2102 forces[:] *= 0
2103 return
2104 dist /= d
2105 df = np.subtract.reduce(forces[self.indices])
2106 f = df.dot(dist)
2107 con_saved = atoms.constraints
2108 try:
2109 con = [con for con in con_saved
2110 if not isinstance(con, MirrorForce)]
2111 atoms.set_constraint(con)
2112 forces_copy = atoms.get_forces()
2113 finally:
2114 atoms.set_constraint(con_saved)
2115 df1 = -1 / 2. * f * dist
2116 forces_copy[self.indices] += (df1, -df1)
2117 # Check if forces would be converged if the bond with mirrored forces
2118 # would also be fixed
2119 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2120 factor = 1.
2121 else:
2122 factor = 0.
2123 df1 = -(1 + factor) / 2. * f * dist
2124 forces[self.indices] += (df1, -df1)
2126 def index_shuffle(self, atoms, ind):
2127 """Shuffle the indices of the two atoms in this constraint
2129 """
2130 newa = [-1, -1] # Signal error
2131 for new, old in slice2enlist(ind, len(atoms)):
2132 for i, a in enumerate(self.indices):
2133 if old == a:
2134 newa[i] = new
2135 if newa[0] == -1 or newa[1] == -1:
2136 raise IndexError('Constraint not part of slice')
2137 self.indices = newa
2139 def __repr__(self):
2140 return 'MirrorForce(%d, %d, %f, %f, %f)' % (
2141 self.indices[0], self.indices[1], self.max_dist, self.min_dist,
2142 self.fmax)
2144 def todict(self):
2145 return {'name': 'MirrorForce',
2146 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2147 'max_dist': self.max_dist,
2148 'min_dist': self.min_dist, 'fmax': self.fmax}}
2151class MirrorTorque(FixConstraint):
2152 """Constraint object for mirroring the torque acting on a dihedral
2153 angle defined by four atoms.
2155 This class is designed to find a transition state with the help of a
2156 single optimization. It can be used if the transition state belongs to a
2157 cis-trans-isomerization with a change of dihedral angle. First the given
2158 dihedral angle will be fixed until all other degrees of freedom are
2159 optimized, then the torque acting on the dihedral angle will be mirrored
2160 to find the transition state. Transition states in
2161 dependence of the force can be obtained by stretching the molecule and
2162 fixing its total length with *FixBondLength* or by using *ExternalForce*
2163 during the optimization with *MirrorTorque*.
2165 This constraint can be used to find
2166 transition states of cis-trans-isomerization.
2168 a1 a4
2169 | |
2170 a2 __ a3
2172 Parameters
2173 ----------
2174 a1: int
2175 First atom index.
2176 a2: int
2177 Second atom index.
2178 a3: int
2179 Third atom index.
2180 a4: int
2181 Fourth atom index.
2182 max_angle: float
2183 Upper limit of the dihedral angle interval where the transition state
2184 can be found.
2185 min_angle: float
2186 Lower limit of the dihedral angle interval where the transition state
2187 can be found.
2188 fmax: float
2189 Maximum force used for the optimization.
2191 Notes
2192 -----
2193 You can combine this constraint for example with FixBondLength but make
2194 sure that *MirrorTorque* comes first in the list if there are overlaps
2195 between atom1-4 and atom5-6:
2197 >>> from ase.build import bulk
2199 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2200 >>> atom1, atom2, atom3, atom4, atom5, atom6 = atoms[:6]
2201 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)
2202 >>> con2 = FixBondLength(atom5, atom6)
2203 >>> atoms.set_constraint([con1, con2])
2205 """
2207 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,
2208 fmax=0.1):
2209 self.indices = [a1, a2, a3, a4]
2210 self.min_angle = min_angle
2211 self.max_angle = max_angle
2212 self.fmax = fmax
2214 def adjust_positions(self, atoms, new):
2215 pass
2217 def adjust_forces(self, atoms, forces):
2218 angle = atoms.get_dihedral(self.indices[0], self.indices[1],
2219 self.indices[2], self.indices[3])
2220 angle *= np.pi / 180.
2221 if (angle < self.min_angle) or (angle > self.max_angle):
2222 # Stop structure optimization
2223 forces[:] *= 0
2224 return
2225 p = atoms.positions[self.indices]
2226 f = forces[self.indices]
2228 f0 = (f[1] + f[2]) / 2.
2229 ff = f - f0
2230 p0 = (p[2] + p[1]) / 2.
2231 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)
2232 fff = ff - np.cross(m0, p - p0)
2233 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \
2234 (p[1] - p0).dot(p[1] - p0)
2235 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \
2236 (p[2] - p0).dot(p[2] - p0)
2237 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \
2238 np.linalg.norm(p[1] - p0)
2239 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \
2240 np.linalg.norm(p[2] - p0)
2241 omegap = omegap1 + omegap2
2242 con_saved = atoms.constraints
2243 try:
2244 con = [con for con in con_saved
2245 if not isinstance(con, MirrorTorque)]
2246 atoms.set_constraint(con)
2247 forces_copy = atoms.get_forces()
2248 finally:
2249 atoms.set_constraint(con_saved)
2250 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \
2251 np.linalg.norm(p[1] - p0)
2252 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \
2253 np.linalg.norm(p[2] - p0)
2254 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2255 # Check if forces would be converged if the dihedral angle with
2256 # mirrored torque would also be fixed
2257 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2258 factor = 1.
2259 else:
2260 factor = 0.
2261 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \
2262 np.linalg.norm(p[1] - p0)
2263 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \
2264 np.linalg.norm(p[2] - p0)
2265 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2267 def index_shuffle(self, atoms, ind):
2268 # See docstring of superclass
2269 indices = []
2270 for new, old in slice2enlist(ind, len(atoms)):
2271 if old in self.indices:
2272 indices.append(new)
2273 if len(indices) == 0:
2274 raise IndexError('All indices in MirrorTorque not part of slice')
2275 self.indices = np.asarray(indices, int)
2277 def __repr__(self):
2278 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (
2279 self.indices[0], self.indices[1], self.indices[2],
2280 self.indices[3], self.max_angle, self.min_angle, self.fmax)
2282 def todict(self):
2283 return {'name': 'MirrorTorque',
2284 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2285 'a3': self.indices[2], 'a4': self.indices[3],
2286 'max_angle': self.max_angle,
2287 'min_angle': self.min_angle, 'fmax': self.fmax}}
2290class FixSymmetry(FixConstraint):
2291 """
2292 Constraint to preserve spacegroup symmetry during optimisation.
2294 Requires spglib package to be available.
2295 """
2297 def __init__(self, atoms, symprec=0.01, adjust_positions=True,
2298 adjust_cell=True, verbose=False):
2299 self.atoms = atoms.copy()
2300 self.symprec = symprec
2301 self.verbose = verbose
2302 refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry
2303 sym = prep_symmetry(atoms, symprec, self.verbose)
2304 self.rotations, self.translations, self.symm_map = sym
2305 self.do_adjust_positions = adjust_positions
2306 self.do_adjust_cell = adjust_cell
2308 def adjust_cell(self, atoms, cell):
2309 if not self.do_adjust_cell:
2310 return
2311 # stress should definitely be symmetrized as a rank 2 tensor
2312 # UnitCellFilter uses deformation gradient as cell DOF with steps
2313 # dF = stress.F^-T quantity that should be symmetrized is therefore dF .
2314 # F^T assume prev F = I, so just symmetrize dF
2315 cur_cell = atoms.get_cell()
2316 cur_cell_inv = atoms.cell.reciprocal().T
2318 # F defined such that cell = cur_cell . F^T
2319 # assume prev F = I, so dF = F - I
2320 delta_deform_grad = np.dot(cur_cell_inv, cell).T - np.eye(3)
2322 # symmetrization doesn't work properly with large steps, since
2323 # it depends on current cell, and cell is being changed by deformation
2324 # gradient
2325 max_delta_deform_grad = np.max(np.abs(delta_deform_grad))
2326 if max_delta_deform_grad > 0.25:
2327 raise RuntimeError('FixSymmetry adjust_cell does not work properly'
2328 ' with large deformation gradient step {} > 0.25'
2329 .format(max_delta_deform_grad))
2330 elif max_delta_deform_grad > 0.15:
2331 warn('FixSymmetry adjust_cell may be ill behaved with large '
2332 'deformation gradient step {}'.format(max_delta_deform_grad))
2334 symmetrized_delta_deform_grad = symmetrize_rank2(cur_cell, cur_cell_inv,
2335 delta_deform_grad,
2336 self.rotations)
2337 cell[:] = np.dot(cur_cell,
2338 (symmetrized_delta_deform_grad + np.eye(3)).T)
2340 def adjust_positions(self, atoms, new):
2341 if not self.do_adjust_positions:
2342 return
2343 # symmetrize changes in position as rank 1 tensors
2344 step = new - atoms.positions
2345 symmetrized_step = symmetrize_rank1(atoms.get_cell(),
2346 atoms.cell.reciprocal().T, step,
2347 self.rotations, self.translations,
2348 self.symm_map)
2349 new[:] = atoms.positions + symmetrized_step
2351 def adjust_forces(self, atoms, forces):
2352 # symmetrize forces as rank 1 tensors
2353 # print('adjusting forces')
2354 forces[:] = symmetrize_rank1(atoms.get_cell(),
2355 atoms.cell.reciprocal().T, forces,
2356 self.rotations, self.translations,
2357 self.symm_map)
2359 def adjust_stress(self, atoms, stress):
2360 # symmetrize stress as rank 2 tensor
2361 raw_stress = voigt_6_to_full_3x3_stress(stress)
2362 symmetrized_stress = symmetrize_rank2(atoms.get_cell(),
2363 atoms.cell.reciprocal().T,
2364 raw_stress, self.rotations)
2365 stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress)
2367 def index_shuffle(self, atoms, ind):
2368 if len(atoms) != len(ind) or len(set(ind)) != len(ind):
2369 raise RuntimeError("FixSymmetry can only accomodate atom"
2370 " permutions, and len(Atoms) == {} "
2371 "!= len(ind) == {} or ind has duplicates"
2372 .format(len(atoms), len(ind)))
2374 ind_reversed = np.zeros((len(ind)), dtype=int)
2375 ind_reversed[ind] = range(len(ind))
2376 new_symm_map = []
2377 for sm in self.symm_map:
2378 new_sm = np.array([-1] * len(atoms))
2379 for at_i in range(len(ind)):
2380 new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]]
2381 new_symm_map.append(new_sm)
2383 self.symm_map = new_symm_map
2385 def todict(self):
2386 return {
2387 'name': 'FixSymmetry',
2388 'kwargs': {
2389 'atoms': self.atoms,
2390 'symprec': self.symprec,
2391 'adjust_positions': self.do_adjust_positions,
2392 'adjust_cell': self.do_adjust_cell,
2393 'verbose': self.verbose,
2394 },
2395 }
2398class Filter(FilterOld):
2399 @deprecated('Import Filter from ase.filters')
2400 def __init__(self, *args, **kwargs):
2401 """
2402 .. deprecated:: 3.23.0
2403 Import ``Filter`` from :mod:`ase.filters`
2404 """
2405 super().__init__(*args, **kwargs)
2408class StrainFilter(StrainFilterOld):
2409 @deprecated('Import StrainFilter from ase.filters')
2410 def __init__(self, *args, **kwargs):
2411 """
2412 .. deprecated:: 3.23.0
2413 Import ``StrainFilter`` from :mod:`ase.filters`
2414 """
2415 super().__init__(*args, **kwargs)
2418class UnitCellFilter(UnitCellFilterOld):
2419 @deprecated('Import UnitCellFilter from ase.filters')
2420 def __init__(self, *args, **kwargs):
2421 """
2422 .. deprecated:: 3.23.0
2423 Import ``UnitCellFilter`` from :mod:`ase.filters`
2424 """
2425 super().__init__(*args, **kwargs)
2428class ExpCellFilter(ExpCellFilterOld):
2429 @deprecated('Import ExpCellFilter from ase.filters')
2430 def __init__(self, *args, **kwargs):
2431 """
2432 .. deprecated:: 3.23.0
2433 Import ``ExpCellFilter`` from :mod:`ase.filters`
2434 or use :class:`~ase.filters.FrechetCellFilter` for better
2435 convergence w.r.t. cell variables
2436 """
2437 super().__init__(*args, **kwargs)