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