Coverage for /builds/hweiske/ase/ase/constraints.py: 86.95%
1218 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +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."""
87 def adjust_positions(self, atoms: Atoms, new):
88 """Adjust positions."""
90 def adjust_momenta(self, atoms: Atoms, momenta):
91 """Adjust momenta."""
92 # The default is in identical manner to forces.
93 # TODO: The default is however not always reasonable.
94 self.adjust_forces(atoms, momenta)
96 def adjust_forces(self, atoms: Atoms, forces):
97 """Adjust forces."""
99 def copy(self):
100 """Copy constraint."""
101 return dict2constraint(self.todict().copy())
103 def todict(self):
104 """Convert constraint to dictionary."""
107class IndexedConstraint(FixConstraint):
108 def __init__(self, indices=None, mask=None):
109 """Constrain chosen atoms.
111 Parameters
112 ----------
113 indices : sequence of int
114 Indices for those atoms that should be constrained.
115 mask : sequence of bool
116 One boolean per atom indicating if the atom should be
117 constrained or not.
118 """
120 if mask is not None:
121 if indices is not None:
122 raise ValueError('Use only one of "indices" and "mask".')
123 indices = mask
124 indices = np.atleast_1d(indices)
125 if np.ndim(indices) > 1:
126 raise ValueError('indices has wrong amount of dimensions. '
127 f'Got {np.ndim(indices)}, expected ndim <= 1')
129 if indices.dtype == bool:
130 indices = np.arange(len(indices))[indices]
131 elif len(indices) == 0:
132 indices = np.empty(0, dtype=int)
133 elif not np.issubdtype(indices.dtype, np.integer):
134 raise ValueError('Indices must be integers or boolean mask, '
135 f'not dtype={indices.dtype}')
137 if len(set(indices)) < len(indices):
138 raise ValueError(
139 'The indices array contains duplicates. '
140 'Perhaps you want to specify a mask instead, but '
141 'forgot the mask= keyword.')
143 self.index = indices
145 def index_shuffle(self, atoms, ind):
146 # See docstring of superclass
147 index = []
149 # Resolve negative indices:
150 actual_indices = set(np.arange(len(atoms))[self.index])
152 for new, old in slice2enlist(ind, len(atoms)):
153 if old in actual_indices:
154 index.append(new)
155 if len(index) == 0:
156 raise IndexError('All indices in FixAtoms not part of slice')
157 self.index = np.asarray(index, int)
158 # XXX make immutable
160 def get_indices(self):
161 return self.index.copy()
163 def repeat(self, m, n):
164 i0 = 0
165 natoms = 0
166 if isinstance(m, int):
167 m = (m, m, m)
168 index_new = []
169 for _ in range(m[2]):
170 for _ in range(m[1]):
171 for _ in range(m[0]):
172 i1 = i0 + n
173 index_new += [i + natoms for i in self.index]
174 i0 = i1
175 natoms += n
176 self.index = np.asarray(index_new, int)
177 # XXX make immutable
178 return self
180 def delete_atoms(self, indices, natoms):
181 """Removes atoms from the index array, if present.
183 Required for removing atoms with existing constraint.
184 """
186 i = np.zeros(natoms, int) - 1
187 new = np.delete(np.arange(natoms), indices)
188 i[new] = np.arange(len(new))
189 index = i[self.index]
190 self.index = index[index >= 0]
191 # XXX make immutable
192 if len(self.index) == 0:
193 return None
194 return self
197class FixAtoms(IndexedConstraint):
198 """Fix chosen atoms.
200 Examples
201 --------
202 Fix all Copper atoms:
204 >>> from ase.build import bulk
206 >>> atoms = bulk('Cu', 'fcc', a=3.6)
207 >>> mask = (atoms.symbols == 'Cu')
208 >>> c = FixAtoms(mask=mask)
209 >>> atoms.set_constraint(c)
211 Fix all atoms with z-coordinate less than 1.0 Angstrom:
213 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)
214 >>> atoms.set_constraint(c)
215 """
217 def get_removed_dof(self, atoms):
218 return 3 * len(self.index)
220 def adjust_positions(self, atoms, new):
221 new[self.index] = atoms.positions[self.index]
223 def adjust_forces(self, atoms, forces):
224 forces[self.index] = 0.0
226 def __repr__(self):
227 clsname = type(self).__name__
228 indices = ints2string(self.index)
229 return f'{clsname}(indices={indices})'
231 def todict(self):
232 return {'name': 'FixAtoms',
233 'kwargs': {'indices': self.index.tolist()}}
236class FixCom(FixConstraint):
237 """Constraint class for fixing the center of mass."""
239 index = slice(None) # all atoms
241 def get_removed_dof(self, atoms):
242 return 3
244 def adjust_positions(self, atoms, new):
245 masses = atoms.get_masses()[self.index]
246 old_cm = atoms.get_center_of_mass(indices=self.index)
247 new_cm = masses @ new[self.index] / masses.sum()
248 diff = old_cm - new_cm
249 new += diff
251 def adjust_momenta(self, atoms, momenta):
252 """Adjust momenta so that the center-of-mass velocity is zero."""
253 masses = atoms.get_masses()[self.index]
254 velocity_com = momenta[self.index].sum(axis=0) / masses.sum()
255 momenta[self.index] -= masses[:, None] * velocity_com
257 def adjust_forces(self, atoms, forces):
258 # Eqs. (3) and (7) in https://doi.org/10.1021/jp9722824
259 masses = atoms.get_masses()[self.index]
260 lmd = masses @ forces[self.index] / sum(masses**2)
261 forces[self.index] -= masses[:, None] * lmd
263 def todict(self):
264 return {'name': 'FixCom',
265 'kwargs': {}}
268class FixSubsetCom(FixCom, IndexedConstraint):
269 """Constraint class for fixing the center of mass of a subset of atoms."""
271 def __init__(self, indices):
272 super().__init__(indices=indices)
274 def todict(self):
275 return {'name': self.__class__.__name__,
276 'kwargs': {'indices': self.index.tolist()}}
279def ints2string(x, threshold=None):
280 """Convert ndarray of ints to string."""
281 if threshold is None or len(x) <= threshold:
282 return str(x.tolist())
283 return str(x[:threshold].tolist())[:-1] + ', ...]'
286class FixBondLengths(FixConstraint):
287 maxiter = 500
289 def __init__(self, pairs, tolerance=1e-13,
290 bondlengths=None, iterations=None):
291 """iterations:
292 Ignored"""
293 self.pairs = np.asarray(pairs)
294 self.tolerance = tolerance
295 self.bondlengths = bondlengths
297 def get_removed_dof(self, atoms):
298 return len(self.pairs)
300 def adjust_positions(self, atoms, new):
301 old = atoms.positions
302 masses = atoms.get_masses()
304 if self.bondlengths is None:
305 self.bondlengths = self.initialize_bond_lengths(atoms)
307 for i in range(self.maxiter):
308 converged = True
309 for j, ab in enumerate(self.pairs):
310 a = ab[0]
311 b = ab[1]
312 cd = self.bondlengths[j]
313 r0 = old[a] - old[b]
314 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
315 d1 = new[a] - new[b] - r0 + d0
316 m = 1 / (1 / masses[a] + 1 / masses[b])
317 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)
318 if abs(x) > self.tolerance:
319 new[a] += x * m / masses[a] * d0
320 new[b] -= x * m / masses[b] * d0
321 converged = False
322 if converged:
323 break
324 else:
325 raise RuntimeError('Did not converge')
327 def adjust_momenta(self, atoms, p):
328 old = atoms.positions
329 masses = atoms.get_masses()
331 if self.bondlengths is None:
332 self.bondlengths = self.initialize_bond_lengths(atoms)
334 for i in range(self.maxiter):
335 converged = True
336 for j, ab in enumerate(self.pairs):
337 a = ab[0]
338 b = ab[1]
339 cd = self.bondlengths[j]
340 d = old[a] - old[b]
341 d, _ = find_mic(d, atoms.cell, atoms.pbc)
342 dv = p[a] / masses[a] - p[b] / masses[b]
343 m = 1 / (1 / masses[a] + 1 / masses[b])
344 x = -np.dot(dv, d) / cd**2
345 if abs(x) > self.tolerance:
346 p[a] += x * m * d
347 p[b] -= x * m * d
348 converged = False
349 if converged:
350 break
351 else:
352 raise RuntimeError('Did not converge')
354 def adjust_forces(self, atoms, forces):
355 self.constraint_forces = -forces
356 self.adjust_momenta(atoms, forces)
357 self.constraint_forces += forces
359 def initialize_bond_lengths(self, atoms):
360 bondlengths = np.zeros(len(self.pairs))
362 for i, ab in enumerate(self.pairs):
363 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)
365 return bondlengths
367 def get_indices(self):
368 return np.unique(self.pairs.ravel())
370 def todict(self):
371 return {'name': 'FixBondLengths',
372 'kwargs': {'pairs': self.pairs.tolist(),
373 'tolerance': self.tolerance}}
375 def index_shuffle(self, atoms, ind):
376 """Shuffle the indices of the two atoms in this constraint"""
377 map = np.zeros(len(atoms), int)
378 map[ind] = 1
379 n = map.sum()
380 map[:] = -1
381 map[ind] = range(n)
382 pairs = map[self.pairs]
383 self.pairs = pairs[(pairs != -1).all(1)]
384 if len(self.pairs) == 0:
385 raise IndexError('Constraint not part of slice')
388def FixBondLength(a1, a2):
389 """Fix distance between atoms with indices a1 and a2."""
390 return FixBondLengths([(a1, a2)])
393class FixLinearTriatomic(FixConstraint):
394 """Holonomic constraints for rigid linear triatomic molecules."""
396 def __init__(self, triples):
397 """Apply RATTLE-type bond constraints between outer atoms n and m
398 and linear vectorial constraints to the position of central
399 atoms o to fix the geometry of linear triatomic molecules of the
400 type:
402 n--o--m
404 Parameters:
406 triples: list
407 Indices of the atoms forming the linear molecules to constrain
408 as triples. Sequence should be (n, o, m) or (m, o, n).
410 When using these constraints in molecular dynamics or structure
411 optimizations, atomic forces need to be redistributed within a
412 triple. The function redistribute_forces_optimization implements
413 the redistribution of forces for structure optimization, while
414 the function redistribute_forces_md implements the redistribution
415 for molecular dynamics.
417 References:
419 Ciccotti et al. Molecular Physics 47 (1982)
420 :doi:`10.1080/00268978200100942`
421 """
422 self.triples = np.asarray(triples)
423 if self.triples.shape[1] != 3:
424 raise ValueError('"triples" has wrong size')
425 self.bondlengths = None
427 def get_removed_dof(self, atoms):
428 return 4 * len(self.triples)
430 @property
431 def n_ind(self):
432 return self.triples[:, 0]
434 @property
435 def m_ind(self):
436 return self.triples[:, 2]
438 @property
439 def o_ind(self):
440 return self.triples[:, 1]
442 def initialize(self, atoms):
443 masses = atoms.get_masses()
444 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)
446 self.bondlengths = self.initialize_bond_lengths(atoms)
447 self.bondlengths_nm = self.bondlengths.sum(axis=1)
449 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]
450 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +
451 C1[:, 1] ** 2 * self.mass_n * self.mass_o +
452 self.mass_n * self.mass_m)
453 C2 = C1 / C2[:, None]
454 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]
455 C3 = C2 * self.mass_o[:, None] * C3[:, None]
456 C3[:, 1] *= -1
457 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T
458 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)
459 C4 = C1 / C4[:, None]
461 self.C1 = C1
462 self.C2 = C2
463 self.C3 = C3
464 self.C4 = C4
466 def adjust_positions(self, atoms, new):
467 old = atoms.positions
468 new_n, new_m, new_o = self.get_slices(new)
470 if self.bondlengths is None:
471 self.initialize(atoms)
473 r0 = old[self.n_ind] - old[self.m_ind]
474 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
475 d1 = new_n - new_m - r0 + d0
476 a = np.einsum('ij,ij->i', d0, d0)
477 b = np.einsum('ij,ij->i', d1, d0)
478 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2
479 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))
480 g = g[:, None] * self.C3
481 new_n -= g[:, 0, None] * d0
482 new_m += g[:, 1, None] * d0
483 if np.allclose(d0, r0):
484 new_o = (self.C1[:, 0, None] * new_n
485 + self.C1[:, 1, None] * new_m)
486 else:
487 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)
488 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)
489 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2
490 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)
492 self.set_slices(new_n, new_m, new_o, new)
494 def adjust_momenta(self, atoms, p):
495 old = atoms.positions
496 p_n, p_m, p_o = self.get_slices(p)
498 if self.bondlengths is None:
499 self.initialize(atoms)
501 mass_nn = self.mass_n[:, None]
502 mass_mm = self.mass_m[:, None]
503 mass_oo = self.mass_o[:, None]
505 d = old[self.n_ind] - old[self.m_ind]
506 d, _ = find_mic(d, atoms.cell, atoms.pbc)
507 dv = p_n / mass_nn - p_m / mass_mm
508 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2
509 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]
510 p_n -= k[:, 0, None] * mass_nn * d
511 p_m += k[:, 1, None] * mass_mm * d
512 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +
513 self.C1[:, 1, None] * p_m / mass_mm))
515 self.set_slices(p_n, p_m, p_o, p)
517 def adjust_forces(self, atoms, forces):
519 if self.bondlengths is None:
520 self.initialize(atoms)
522 A = self.C4 * np.diff(self.C1)
523 A[:, 0] *= -1
524 A -= 1
525 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]
526 A /= (A.sum(axis=1))[:, None]
528 self.constraint_forces = -forces
529 old = atoms.positions
531 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)
533 d = old[self.n_ind] - old[self.m_ind]
534 d, _ = find_mic(d, atoms.cell, atoms.pbc)
535 df = fr_n - fr_m
536 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2
537 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]
538 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]
539 forces[self.o_ind] = fr_o + k[:, None] * d * B
541 self.constraint_forces += forces
543 def redistribute_forces_optimization(self, forces):
544 """Redistribute forces within a triple when performing structure
545 optimizations.
547 The redistributed forces needs to be further adjusted using the
548 appropriate Lagrange multipliers as implemented in adjust_forces."""
549 forces_n, forces_m, forces_o = self.get_slices(forces)
550 C1_1 = self.C1[:, 0, None]
551 C1_2 = self.C1[:, 1, None]
552 C4_1 = self.C4[:, 0, None]
553 C4_2 = self.C4[:, 1, None]
555 fr_n = ((1 - C4_1 * C1_1) * forces_n -
556 C4_1 * (C1_2 * forces_m - forces_o))
557 fr_m = ((1 - C4_2 * C1_2) * forces_m -
558 C4_2 * (C1_1 * forces_n - forces_o))
559 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +
560 C4_1 * forces_n + C4_2 * forces_m)
562 return fr_n, fr_m, fr_o
564 def redistribute_forces_md(self, atoms, forces, rand=False):
565 """Redistribute forces within a triple when performing molecular
566 dynamics.
568 When rand=True, use the equations for random force terms, as
569 used e.g. by Langevin dynamics, otherwise apply the standard
570 equations for deterministic forces (see Ciccotti et al. Molecular
571 Physics 47 (1982))."""
572 if self.bondlengths is None:
573 self.initialize(atoms)
574 forces_n, forces_m, forces_o = self.get_slices(forces)
575 C1_1 = self.C1[:, 0, None]
576 C1_2 = self.C1[:, 1, None]
577 C2_1 = self.C2[:, 0, None]
578 C2_2 = self.C2[:, 1, None]
579 mass_nn = self.mass_n[:, None]
580 mass_mm = self.mass_m[:, None]
581 mass_oo = self.mass_o[:, None]
582 if rand:
583 mr1 = (mass_mm / mass_nn) ** 0.5
584 mr2 = (mass_oo / mass_nn) ** 0.5
585 mr3 = (mass_nn / mass_mm) ** 0.5
586 mr4 = (mass_oo / mass_mm) ** 0.5
587 else:
588 mr1 = 1.0
589 mr2 = 1.0
590 mr3 = 1.0
591 mr4 = 1.0
593 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -
594 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -
595 mr2 * mass_mm * mass_nn * forces_o))
597 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -
598 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -
599 mr4 * mass_mm * mass_nn * forces_o))
601 self.set_slices(fr_n, fr_m, 0.0, forces)
603 def get_slices(self, a):
604 a_n = a[self.n_ind]
605 a_m = a[self.m_ind]
606 a_o = a[self.o_ind]
608 return a_n, a_m, a_o
610 def set_slices(self, a_n, a_m, a_o, a):
611 a[self.n_ind] = a_n
612 a[self.m_ind] = a_m
613 a[self.o_ind] = a_o
615 def initialize_bond_lengths(self, atoms):
616 bondlengths = np.zeros((len(self.triples), 2))
618 for i in range(len(self.triples)):
619 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],
620 self.o_ind[i], mic=True)
621 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],
622 self.m_ind[i], mic=True)
624 return bondlengths
626 def get_indices(self):
627 return np.unique(self.triples.ravel())
629 def todict(self):
630 return {'name': 'FixLinearTriatomic',
631 'kwargs': {'triples': self.triples.tolist()}}
633 def index_shuffle(self, atoms, ind):
634 """Shuffle the indices of the three atoms in this constraint"""
635 map = np.zeros(len(atoms), int)
636 map[ind] = 1
637 n = map.sum()
638 map[:] = -1
639 map[ind] = range(n)
640 triples = map[self.triples]
641 self.triples = triples[(triples != -1).all(1)]
642 if len(self.triples) == 0:
643 raise IndexError('Constraint not part of slice')
646class FixedMode(FixConstraint):
647 """Constrain atoms to move along directions orthogonal to
648 a given mode only. Initialize with a mode, such as one produced by
649 ase.vibrations.Vibrations.get_mode()."""
651 def __init__(self, mode):
652 mode = np.asarray(mode)
653 self.mode = (mode / np.sqrt((mode**2).sum())).reshape(-1)
655 def get_removed_dof(self, atoms):
656 return len(atoms)
658 def adjust_positions(self, atoms, newpositions):
659 newpositions = newpositions.ravel()
660 oldpositions = atoms.positions.ravel()
661 step = newpositions - oldpositions
662 newpositions -= self.mode * np.dot(step, self.mode)
664 def adjust_forces(self, atoms, forces):
665 forces = forces.ravel()
666 forces -= self.mode * np.dot(forces, self.mode)
668 def index_shuffle(self, atoms, ind):
669 eps = 1e-12
670 mode = self.mode.reshape(-1, 3)
671 excluded = np.ones(len(mode), dtype=bool)
672 excluded[ind] = False
673 if (abs(mode[excluded]) > eps).any():
674 raise IndexError('All nonzero parts of mode not in slice')
675 self.mode = mode[ind].ravel()
677 def get_indices(self):
678 # This function will never properly work because it works on all
679 # atoms and it has no idea how to tell how many atoms it is
680 # attached to. If it is being used, surely the user knows
681 # everything is being constrained.
682 return []
684 def todict(self):
685 return {'name': 'FixedMode',
686 'kwargs': {'mode': self.mode.tolist()}}
688 def __repr__(self):
689 return f'FixedMode({self.mode.tolist()})'
692def _normalize(direction):
693 if np.shape(direction) != (3,):
694 raise ValueError("len(direction) is {len(direction)}. Has to be 3")
696 direction = np.asarray(direction) / np.linalg.norm(direction)
697 return direction
700class FixedPlane(IndexedConstraint):
701 """
702 Constraint object for fixing chosen atoms to only move in a plane.
704 The plane is defined by its normal vector *direction*
705 """
707 def __init__(self, indices, direction):
708 """Constrain chosen atoms.
710 Parameters
711 ----------
712 indices : int or list of int
713 Index or indices for atoms that should be constrained
714 direction : list of 3 int
715 Direction of the normal vector
717 Examples
718 --------
719 Fix all Copper atoms to only move in the yz-plane:
721 >>> from ase.build import bulk
722 >>> from ase.constraints import FixedPlane
724 >>> atoms = bulk('Cu', 'fcc', a=3.6)
725 >>> c = FixedPlane(
726 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
727 ... direction=[1, 0, 0],
728 ... )
729 >>> atoms.set_constraint(c)
731 or constrain a single atom with the index 0 to move in the xy-plane:
733 >>> c = FixedPlane(indices=0, direction=[0, 0, 1])
734 >>> atoms.set_constraint(c)
735 """
736 super().__init__(indices=indices)
737 self.dir = _normalize(direction)
739 def adjust_positions(self, atoms, newpositions):
740 step = newpositions[self.index] - atoms.positions[self.index]
741 newpositions[self.index] -= _projection(step, self.dir)
743 def adjust_forces(self, atoms, forces):
744 forces[self.index] -= _projection(forces[self.index], self.dir)
746 def get_removed_dof(self, atoms):
747 return len(self.index)
749 def todict(self):
750 return {
751 'name': 'FixedPlane',
752 'kwargs': {'indices': self.index.tolist(),
753 'direction': self.dir.tolist()}
754 }
756 def __repr__(self):
757 return f'FixedPlane(indices={self.index}, {self.dir.tolist()})'
760def _projection(vectors, direction):
761 dotprods = vectors @ direction
762 projection = direction[None, :] * dotprods[:, None]
763 return projection
766class FixedLine(IndexedConstraint):
767 """
768 Constrain an atom index or a list of atom indices to move on a line only.
770 The line is defined by its vector *direction*
771 """
773 def __init__(self, indices, direction):
774 """Constrain chosen atoms.
776 Parameters
777 ----------
778 indices : int or list of int
779 Index or indices for atoms that should be constrained
780 direction : list of 3 int
781 Direction of the vector defining the line
783 Examples
784 --------
785 Fix all Copper atoms to only move in the x-direction:
787 >>> from ase.constraints import FixedLine
788 >>> c = FixedLine(
789 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
790 ... direction=[1, 0, 0],
791 ... )
792 >>> atoms.set_constraint(c)
794 or constrain a single atom with the index 0 to move in the z-direction:
796 >>> c = FixedLine(indices=0, direction=[0, 0, 1])
797 >>> atoms.set_constraint(c)
798 """
799 super().__init__(indices)
800 self.dir = _normalize(direction)
802 def adjust_positions(self, atoms, newpositions):
803 step = newpositions[self.index] - atoms.positions[self.index]
804 projection = _projection(step, self.dir)
805 newpositions[self.index] = atoms.positions[self.index] + projection
807 def adjust_forces(self, atoms, forces):
808 forces[self.index] = _projection(forces[self.index], self.dir)
810 def get_removed_dof(self, atoms):
811 return 2 * len(self.index)
813 def __repr__(self):
814 return f'FixedLine(indices={self.index}, {self.dir.tolist()})'
816 def todict(self):
817 return {
818 'name': 'FixedLine',
819 'kwargs': {'indices': self.index.tolist(),
820 'direction': self.dir.tolist()}
821 }
824class FixCartesian(IndexedConstraint):
825 """Fix atoms in the directions of the cartesian coordinates.
827 Parameters
828 ----------
829 a : Sequence[int]
830 Indices of atoms to be fixed.
831 mask : tuple[bool, bool, bool], default: (True, True, True)
832 Cartesian directions to be fixed. (False: unfixed, True: fixed)
833 """
835 def __init__(self, a, mask=(True, True, True)):
836 super().__init__(indices=a)
837 self.mask = np.asarray(mask, bool)
839 def get_removed_dof(self, atoms: Atoms):
840 return self.mask.sum() * len(self.index)
842 def adjust_positions(self, atoms: Atoms, new):
843 new[self.index] = np.where(
844 self.mask[None, :],
845 atoms.positions[self.index],
846 new[self.index],
847 )
849 def adjust_forces(self, atoms: Atoms, forces):
850 forces[self.index] *= ~self.mask[None, :]
852 def todict(self):
853 return {'name': 'FixCartesian',
854 'kwargs': {'a': self.index.tolist(),
855 'mask': self.mask.tolist()}}
857 def __repr__(self):
858 name = type(self).__name__
859 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'
862class FixScaled(IndexedConstraint):
863 """Fix atoms in the directions of the unit vectors.
865 Parameters
866 ----------
867 a : Sequence[int]
868 Indices of atoms to be fixed.
869 mask : tuple[bool, bool, bool], default: (True, True, True)
870 Cell directions to be fixed. (False: unfixed, True: fixed)
871 """
873 def __init__(self, a, mask=(True, True, True), cell=None):
874 # XXX The unused cell keyword is there for compatibility
875 # with old trajectory files.
876 super().__init__(indices=a)
877 self.mask = np.asarray(mask, bool)
879 def get_removed_dof(self, atoms: Atoms):
880 return self.mask.sum() * len(self.index)
882 def adjust_positions(self, atoms: Atoms, new):
883 cell = atoms.cell
884 scaled_old = cell.scaled_positions(atoms.positions[self.index])
885 scaled_new = cell.scaled_positions(new[self.index])
886 scaled_new[:, self.mask] = scaled_old[:, self.mask]
887 new[self.index] = cell.cartesian_positions(scaled_new)
889 def adjust_forces(self, atoms: Atoms, forces):
890 # Forces are covariant to the coordinate transformation,
891 # use the inverse transformations
892 cell = atoms.cell
893 scaled_forces = cell.cartesian_positions(forces[self.index])
894 scaled_forces *= -(self.mask - 1)
895 forces[self.index] = cell.scaled_positions(scaled_forces)
897 def todict(self):
898 return {'name': 'FixScaled',
899 'kwargs': {'a': self.index.tolist(),
900 'mask': self.mask.tolist()}}
902 def __repr__(self):
903 name = type(self).__name__
904 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'
907# TODO: Better interface might be to use dictionaries in place of very
908# nested lists/tuples
909class FixInternals(FixConstraint):
910 """Constraint object for fixing multiple internal coordinates.
912 Allows fixing bonds, angles, dihedrals as well as linear combinations
913 of bonds (bondcombos).
915 Please provide angular units in degrees using `angles_deg` and
916 `dihedrals_deg`.
917 Fixing planar angles is not supported at the moment.
918 """
920 def __init__(self, bonds=None, angles=None, dihedrals=None,
921 angles_deg=None, dihedrals_deg=None,
922 bondcombos=None,
923 mic=False, epsilon=1.e-7):
924 """
925 A constrained internal coordinate is defined as a nested list:
926 '[value, [atom indices]]'. The constraint is initialized with a list of
927 constrained internal coordinates, i.e. '[[value, [atom indices]], ...]'.
928 If 'value' is None, the current value of the coordinate is constrained.
930 Parameters
931 ----------
932 bonds: nested python list, optional
933 List with targetvalue and atom indices defining the fixed bonds,
934 i.e. [[targetvalue, [index0, index1]], ...]
936 angles_deg: nested python list, optional
937 List with targetvalue and atom indices defining the fixedangles,
938 i.e. [[targetvalue, [index0, index1, index3]], ...]
940 dihedrals_deg: nested python list, optional
941 List with targetvalue and atom indices defining the fixed dihedrals,
942 i.e. [[targetvalue, [index0, index1, index3]], ...]
944 bondcombos: nested python list, optional
945 List with targetvalue, atom indices and linear coefficient defining
946 the fixed linear combination of bonds,
947 i.e. [[targetvalue, [[index0, index1, coefficient_for_bond],
948 [index1, index2, coefficient_for_bond]]], ...]
950 mic: bool, optional, default: False
951 Minimum image convention.
953 epsilon: float, optional, default: 1e-7
954 Convergence criterion.
955 """
956 warn_msg = 'Please specify {} in degrees using the {} argument.'
957 if angles:
958 warn(warn_msg.format('angles', 'angle_deg'), FutureWarning)
959 angles = np.asarray(angles)
960 angles[:, 0] = angles[:, 0] / np.pi * 180
961 angles = angles.tolist()
962 else:
963 angles = angles_deg
964 if dihedrals:
965 warn(warn_msg.format('dihedrals', 'dihedrals_deg'), FutureWarning)
966 dihedrals = np.asarray(dihedrals)
967 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180
968 dihedrals = dihedrals.tolist()
969 else:
970 dihedrals = dihedrals_deg
972 self.bonds = bonds or []
973 self.angles = angles or []
974 self.dihedrals = dihedrals or []
975 self.bondcombos = bondcombos or []
976 self.mic = mic
977 self.epsilon = epsilon
979 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)
980 + len(self.bondcombos))
982 # Initialize these at run-time:
983 self.constraints = []
984 self.initialized = False
986 def get_removed_dof(self, atoms):
987 return self.n
989 def initialize(self, atoms):
990 if self.initialized:
991 return
992 masses = np.repeat(atoms.get_masses(), 3)
993 cell = None
994 pbc = None
995 if self.mic:
996 cell = atoms.cell
997 pbc = atoms.pbc
998 self.constraints = []
999 for data, ConstrClass in [(self.bonds, self.FixBondLengthAlt),
1000 (self.angles, self.FixAngle),
1001 (self.dihedrals, self.FixDihedral),
1002 (self.bondcombos, self.FixBondCombo)]:
1003 for datum in data:
1004 targetvalue = datum[0]
1005 if targetvalue is None: # set to current value
1006 targetvalue = ConstrClass.get_value(atoms, datum[1],
1007 self.mic)
1008 constr = ConstrClass(targetvalue, datum[1], masses, cell, pbc)
1009 self.constraints.append(constr)
1010 self.initialized = True
1012 @staticmethod
1013 def get_bondcombo(atoms, indices, mic=False):
1014 """Convenience function to return the value of the bondcombo coordinate
1015 (linear combination of bond lengths) for the given Atoms object 'atoms'.
1016 Example: Get the current value of the linear combination of two bond
1017 lengths defined as `bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]`."""
1018 c = sum(df[2] * atoms.get_distance(*df[:2], mic=mic) for df in indices)
1019 return c
1021 def get_subconstraint(self, atoms, definition):
1022 """Get pointer to a specific subconstraint.
1023 Identification by its definition via indices (and coefficients)."""
1024 self.initialize(atoms)
1025 for subconstr in self.constraints:
1026 if isinstance(definition[0], Sequence): # Combo constraint
1027 defin = [d + [c] for d, c in zip(subconstr.indices,
1028 subconstr.coefs)]
1029 if defin == definition:
1030 return subconstr
1031 else: # identify primitive constraints by their indices
1032 if subconstr.indices == [definition]:
1033 return subconstr
1034 raise ValueError('Given `definition` not found on Atoms object.')
1036 def shuffle_definitions(self, shuffle_dic, internal_type):
1037 dfns = [] # definitions
1038 for dfn in internal_type: # e.g. for bond in self.bonds
1039 append = True
1040 new_dfn = [dfn[0], list(dfn[1])]
1041 for old in dfn[1]:
1042 if old in shuffle_dic:
1043 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]
1044 else:
1045 append = False
1046 break
1047 if append:
1048 dfns.append(new_dfn)
1049 return dfns
1051 def shuffle_combos(self, shuffle_dic, internal_type):
1052 dfns = [] # definitions
1053 for dfn in internal_type: # i.e. for bondcombo in self.bondcombos
1054 append = True
1055 all_indices = [idx[0:-1] for idx in dfn[1]]
1056 new_dfn = [dfn[0], list(dfn[1])]
1057 for i, indices in enumerate(all_indices):
1058 for old in indices:
1059 if old in shuffle_dic:
1060 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]
1061 else:
1062 append = False
1063 break
1064 if not append:
1065 break
1066 if append:
1067 dfns.append(new_dfn)
1068 return dfns
1070 def index_shuffle(self, atoms, ind):
1071 # See docstring of superclass
1072 self.initialize(atoms)
1073 shuffle_dic = dict(slice2enlist(ind, len(atoms)))
1074 shuffle_dic = {old: new for new, old in shuffle_dic.items()}
1075 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)
1076 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)
1077 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)
1078 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)
1079 self.initialized = False
1080 self.initialize(atoms)
1081 if len(self.constraints) == 0:
1082 raise IndexError('Constraint not part of slice')
1084 def get_indices(self):
1085 cons = []
1086 for dfn in self.bonds + self.dihedrals + self.angles:
1087 cons.extend(dfn[1])
1088 for dfn in self.bondcombos:
1089 for partial_dfn in dfn[1]:
1090 cons.extend(partial_dfn[0:-1]) # last index is the coefficient
1091 return list(set(cons))
1093 def todict(self):
1094 return {'name': 'FixInternals',
1095 'kwargs': {'bonds': self.bonds,
1096 'angles_deg': self.angles,
1097 'dihedrals_deg': self.dihedrals,
1098 'bondcombos': self.bondcombos,
1099 'mic': self.mic,
1100 'epsilon': self.epsilon}}
1102 def adjust_positions(self, atoms, newpos):
1103 self.initialize(atoms)
1104 for constraint in self.constraints:
1105 constraint.setup_jacobian(atoms.positions)
1106 for _ in range(50):
1107 maxerr = 0.0
1108 for constraint in self.constraints:
1109 constraint.adjust_positions(atoms.positions, newpos)
1110 maxerr = max(abs(constraint.sigma), maxerr)
1111 if maxerr < self.epsilon:
1112 return
1113 msg = 'FixInternals.adjust_positions did not converge.'
1114 if any(constr.targetvalue > 175. or constr.targetvalue < 5. for constr
1115 in self.constraints if isinstance(constr, self.FixAngle)):
1116 msg += (' This may be caused by an almost planar angle.'
1117 ' Support for planar angles would require the'
1118 ' implementation of ghost, i.e. dummy, atoms.'
1119 ' See issue #868.')
1120 raise ValueError(msg)
1122 def adjust_forces(self, atoms, forces):
1123 """Project out translations and rotations and all other constraints"""
1124 self.initialize(atoms)
1125 positions = atoms.positions
1126 N = len(forces)
1127 list2_constraints = list(np.zeros((6, N, 3)))
1128 tx, ty, tz, rx, ry, rz = list2_constraints
1130 list_constraints = [r.ravel() for r in list2_constraints]
1132 tx[:, 0] = 1.0
1133 ty[:, 1] = 1.0
1134 tz[:, 2] = 1.0
1135 ff = forces.ravel()
1137 # Calculate the center of mass
1138 center = positions.sum(axis=0) / N
1140 rx[:, 1] = -(positions[:, 2] - center[2])
1141 rx[:, 2] = positions[:, 1] - center[1]
1142 ry[:, 0] = positions[:, 2] - center[2]
1143 ry[:, 2] = -(positions[:, 0] - center[0])
1144 rz[:, 0] = -(positions[:, 1] - center[1])
1145 rz[:, 1] = positions[:, 0] - center[0]
1147 # Normalizing transl., rotat. constraints
1148 for r in list2_constraints:
1149 r /= np.linalg.norm(r.ravel())
1151 # Add all angle, etc. constraint vectors
1152 for constraint in self.constraints:
1153 constraint.setup_jacobian(positions)
1154 constraint.adjust_forces(positions, forces)
1155 list_constraints.insert(0, constraint.jacobian)
1156 # QR DECOMPOSITION - GRAM SCHMIDT
1158 list_constraints = [r.ravel() for r in list_constraints]
1159 aa = np.column_stack(list_constraints)
1160 (aa, bb) = np.linalg.qr(aa)
1161 # Projection
1162 hh = []
1163 for i, constraint in enumerate(self.constraints):
1164 hh.append(aa[:, i] * np.row_stack(aa[:, i]))
1166 txx = aa[:, self.n] * np.row_stack(aa[:, self.n])
1167 tyy = aa[:, self.n + 1] * np.row_stack(aa[:, self.n + 1])
1168 tzz = aa[:, self.n + 2] * np.row_stack(aa[:, self.n + 2])
1169 rxx = aa[:, self.n + 3] * np.row_stack(aa[:, self.n + 3])
1170 ryy = aa[:, self.n + 4] * np.row_stack(aa[:, self.n + 4])
1171 rzz = aa[:, self.n + 5] * np.row_stack(aa[:, self.n + 5])
1172 T = txx + tyy + tzz + rxx + ryy + rzz
1173 for vec in hh:
1174 T += vec
1175 ff = np.dot(T, np.row_stack(ff))
1176 forces[:, :] -= np.dot(T, np.row_stack(ff)).reshape(-1, 3)
1178 def __repr__(self):
1179 constraints = [repr(constr) for constr in self.constraints]
1180 return f'FixInternals(_copy_init={constraints}, epsilon={self.epsilon})'
1182 # Classes for internal use in FixInternals
1183 class FixInternalsBase:
1184 """Base class for subclasses of FixInternals."""
1186 def __init__(self, targetvalue, indices, masses, cell, pbc):
1187 self.targetvalue = targetvalue # constant target value
1188 self.indices = [defin[0:-1] for defin in indices] # indices, defs
1189 self.coefs = np.asarray([defin[-1] for defin in indices])
1190 self.masses = masses
1191 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix
1192 self.sigma = 1. # difference between current and target value
1193 self.projected_force = None # helps optimizers scan along constr.
1194 self.cell = cell
1195 self.pbc = pbc
1197 def finalize_jacobian(self, pos, n_internals, n, derivs):
1198 """Populate jacobian with derivatives for `n_internals` defined
1199 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""
1200 jacobian = np.zeros((n_internals, *pos.shape))
1201 for i, idx in enumerate(self.indices):
1202 for j in range(n):
1203 jacobian[i, idx[j]] = derivs[i, j]
1204 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))
1205 return self.coefs @ jacobian
1207 def finalize_positions(self, newpos):
1208 jacobian = self.jacobian / self.masses
1209 lamda = -self.sigma / (jacobian @ self.get_jacobian(newpos))
1210 dnewpos = lamda * jacobian
1211 newpos += dnewpos.reshape(newpos.shape)
1213 def adjust_forces(self, positions, forces):
1214 self.projected_forces = ((self.jacobian @ forces.ravel())
1215 * self.jacobian)
1216 self.jacobian /= np.linalg.norm(self.jacobian)
1218 class FixBondCombo(FixInternalsBase):
1219 """Constraint subobject for fixing linear combination of bond lengths
1220 within FixInternals.
1222 sum_i( coef_i * bond_length_i ) = constant
1223 """
1225 def get_jacobian(self, pos):
1226 bondvectors = [pos[k] - pos[h] for h, k in self.indices]
1227 derivs = get_distances_derivatives(bondvectors, cell=self.cell,
1228 pbc=self.pbc)
1229 return self.finalize_jacobian(pos, len(bondvectors), 2, derivs)
1231 def setup_jacobian(self, pos):
1232 self.jacobian = self.get_jacobian(pos)
1234 def adjust_positions(self, oldpos, newpos):
1235 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]
1236 (_, ), (dists, ) = conditional_find_mic([bondvectors],
1237 cell=self.cell,
1238 pbc=self.pbc)
1239 value = self.coefs @ dists
1240 self.sigma = value - self.targetvalue
1241 self.finalize_positions(newpos)
1243 @staticmethod
1244 def get_value(atoms, indices, mic):
1245 return FixInternals.get_bondcombo(atoms, indices, mic)
1247 def __repr__(self):
1248 return (f'FixBondCombo({self.targetvalue}, {self.indices}, '
1249 '{self.coefs})')
1251 class FixBondLengthAlt(FixBondCombo):
1252 """Constraint subobject for fixing bond length within FixInternals.
1253 Fix distance between atoms with indices a1, a2."""
1255 def __init__(self, targetvalue, indices, masses, cell, pbc):
1256 if targetvalue <= 0.:
1257 raise ZeroDivisionError('Invalid targetvalue for fixed bond')
1258 indices = [list(indices) + [1.]] # bond definition with coef 1.
1259 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1261 @staticmethod
1262 def get_value(atoms, indices, mic):
1263 return atoms.get_distance(*indices, mic=mic)
1265 def __repr__(self):
1266 return f'FixBondLengthAlt({self.targetvalue}, {self.indices})'
1268 class FixAngle(FixInternalsBase):
1269 """Constraint subobject for fixing an angle within FixInternals.
1271 Convergence is potentially problematic for angles very close to
1272 0 or 180 degrees as there is a singularity in the Cartesian derivative.
1273 Fixing planar angles is therefore not supported at the moment.
1274 """
1276 def __init__(self, targetvalue, indices, masses, cell, pbc):
1277 """Fix atom movement to construct a constant angle."""
1278 if targetvalue <= 0. or targetvalue >= 180.:
1279 raise ZeroDivisionError('Invalid targetvalue for fixed angle')
1280 indices = [list(indices) + [1.]] # angle definition with coef 1.
1281 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1283 def gather_vectors(self, pos):
1284 v0 = [pos[h] - pos[k] for h, k, l in self.indices]
1285 v1 = [pos[l] - pos[k] for h, k, l in self.indices]
1286 return v0, v1
1288 def get_jacobian(self, pos):
1289 v0, v1 = self.gather_vectors(pos)
1290 derivs = get_angles_derivatives(v0, v1, cell=self.cell,
1291 pbc=self.pbc)
1292 return self.finalize_jacobian(pos, len(v0), 3, derivs)
1294 def setup_jacobian(self, pos):
1295 self.jacobian = self.get_jacobian(pos)
1297 def adjust_positions(self, oldpos, newpos):
1298 v0, v1 = self.gather_vectors(newpos)
1299 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)
1300 self.sigma = value - self.targetvalue
1301 self.finalize_positions(newpos)
1303 @staticmethod
1304 def get_value(atoms, indices, mic):
1305 return atoms.get_angle(*indices, mic=mic)
1307 def __repr__(self):
1308 return f'FixAngle({self.targetvalue}, {self.indices})'
1310 class FixDihedral(FixInternalsBase):
1311 """Constraint subobject for fixing a dihedral angle within FixInternals.
1313 A dihedral becomes undefined when at least one of the inner two angles
1314 becomes planar. Make sure to avoid this situation.
1315 """
1317 def __init__(self, targetvalue, indices, masses, cell, pbc):
1318 indices = [list(indices) + [1.]] # dihedral def. with coef 1.
1319 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1321 def gather_vectors(self, pos):
1322 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]
1323 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]
1324 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]
1325 return v0, v1, v2
1327 def get_jacobian(self, pos):
1328 v0, v1, v2 = self.gather_vectors(pos)
1329 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,
1330 pbc=self.pbc)
1331 return self.finalize_jacobian(pos, len(v0), 4, derivs)
1333 def setup_jacobian(self, pos):
1334 self.jacobian = self.get_jacobian(pos)
1336 def adjust_positions(self, oldpos, newpos):
1337 v0, v1, v2 = self.gather_vectors(newpos)
1338 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1339 # apply minimum dihedral difference 'convention': (diff <= 180)
1340 self.sigma = (value - self.targetvalue + 180) % 360 - 180
1341 self.finalize_positions(newpos)
1343 @staticmethod
1344 def get_value(atoms, indices, mic):
1345 return atoms.get_dihedral(*indices, mic=mic)
1347 def __repr__(self):
1348 return f'FixDihedral({self.targetvalue}, {self.indices})'
1351class FixParametricRelations(FixConstraint):
1353 def __init__(
1354 self,
1355 indices,
1356 Jacobian,
1357 const_shift,
1358 params=None,
1359 eps=1e-12,
1360 use_cell=False,
1361 ):
1362 """Constrains the degrees of freedom to act in a reduced parameter
1363 space defined by the Jacobian
1365 These constraints are based off the work in:
1366 https://arxiv.org/abs/1908.01610
1368 The constraints linearly maps the full 3N degrees of freedom,
1369 where N is number of active lattice vectors/atoms onto a
1370 reduced subset of M free parameters, where M <= 3*N. The
1371 Jacobian matrix and constant shift vector map the full set of
1372 degrees of freedom onto the reduced parameter space.
1374 Currently the constraint is set up to handle either atomic
1375 positions or lattice vectors at one time, but not both. To do
1376 both simply add a two constraints for each set. This is done
1377 to keep the mathematics behind the operations separate.
1379 It would be possible to extend these constraints to allow
1380 non-linear transformations if functionality to update the
1381 Jacobian at each position update was included. This would
1382 require passing an update function evaluate it every time
1383 adjust_positions is callled. This is currently NOT supported,
1384 and there are no plans to implement it in the future.
1386 Args:
1387 indices (list of int): indices of the constrained atoms
1388 (if not None or empty then cell_indices must be None or Empty)
1389 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))):
1390 The Jacobian describing
1391 the parameter space transformation
1392 const_shift (np.ndarray(shape=(3*len(indices)))):
1393 A vector describing the constant term
1394 in the transformation not accounted for in the Jacobian
1395 params (list of str):
1396 parameters used in the parametric representation
1397 if None a list is generated based on the shape of the Jacobian
1398 eps (float): a small number to compare the similarity of
1399 numbers and set the precision used
1400 to generate the constraint expressions
1401 use_cell (bool): if True then act on the cell object
1403 """
1404 self.indices = np.array(indices)
1405 self.Jacobian = np.array(Jacobian)
1406 self.const_shift = np.array(const_shift)
1408 assert self.const_shift.shape[0] == 3 * len(self.indices)
1409 assert self.Jacobian.shape[0] == 3 * len(self.indices)
1411 self.eps = eps
1412 self.use_cell = use_cell
1414 if params is None:
1415 params = []
1416 if self.Jacobian.shape[1] > 0:
1417 int_fmt_str = "{:0" + \
1418 str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"
1419 for param_ind in range(self.Jacobian.shape[1]):
1420 params.append("param_" + int_fmt_str.format(param_ind))
1421 else:
1422 assert len(params) == self.Jacobian.shape[-1]
1424 self.params = params
1426 self.Jacobian_inv = np.linalg.inv(
1427 self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T
1429 @classmethod
1430 def from_expressions(cls, indices, params, expressions,
1431 eps=1e-12, use_cell=False):
1432 """Converts the expressions into a Jacobian Matrix/const_shift
1433 vector and constructs a FixParametricRelations constraint
1435 The expressions must be a list like object of size 3*N and
1436 elements must be ordered as:
1437 [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],
1438 where i, j, and k are the first, second and third
1439 component of the atomic position/lattice
1440 vector. Currently only linear operations are allowed to be
1441 included in the expressions so
1442 only terms like:
1443 - const * param_0
1444 - sqrt[const] * param_1
1445 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M
1446 where const is any real number and param_0, param_1, ..., param_M are
1447 the parameters passed in
1448 params, are allowed.
1450 For example, fractional atomic position constraints for wurtzite are:
1451 params = ["z1", "z2"]
1452 expressions = [
1453 "1.0/3.0", "2.0/3.0", "z1",
1454 "2.0/3.0", "1.0/3.0", "0.5 + z1",
1455 "1.0/3.0", "2.0/3.0", "z2",
1456 "2.0/3.0", "1.0/3.0", "0.5 + z2",
1457 ]
1459 For diamond are:
1460 params = []
1461 expressions = [
1462 "0.0", "0.0", "0.0",
1463 "0.25", "0.25", "0.25",
1464 ],
1466 and for stannite are
1467 params=["x4", "z4"]
1468 expressions = [
1469 "0.0", "0.0", "0.0",
1470 "0.0", "0.5", "0.5",
1471 "0.75", "0.25", "0.5",
1472 "0.25", "0.75", "0.5",
1473 "x4 + z4", "x4 + z4", "2*x4",
1474 "x4 - z4", "x4 - z4", "-2*x4",
1475 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",
1476 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",
1477 ]
1479 Args:
1480 indices (list of int): indices of the constrained atoms
1481 (if not None or empty then cell_indices must be None or Empty)
1482 params (list of str): parameters used in the
1483 parametric representation
1484 expressions (list of str): expressions used to convert from the
1485 parametric to the real space representation
1486 eps (float): a small number to compare the similarity of
1487 numbers and set the precision used
1488 to generate the constraint expressions
1489 use_cell (bool): if True then act on the cell object
1491 Returns:
1492 cls(
1493 indices,
1494 Jacobian generated from expressions,
1495 const_shift generated from expressions,
1496 params,
1497 eps-12,
1498 use_cell,
1499 )
1500 """
1501 Jacobian = np.zeros((3 * len(indices), len(params)))
1502 const_shift = np.zeros(3 * len(indices))
1504 for expr_ind, expression in enumerate(expressions):
1505 expression = expression.strip()
1507 # Convert subtraction to addition
1508 expression = expression.replace("-", "+(-1.0)*")
1509 if expression[0] == "+":
1510 expression = expression[1:]
1511 elif expression[:2] == "(+":
1512 expression = "(" + expression[2:]
1514 # Explicitly add leading zeros so when replacing param_1 with 0.0
1515 # param_11 does not become 0.01
1516 int_fmt_str = "{:0" + \
1517 str(int(np.ceil(np.log10(len(params) + 1)))) + "d}"
1519 param_dct = {}
1520 param_map = {}
1522 # Construct a standardized param template for A/B filling
1523 for param_ind, param in enumerate(params):
1524 param_str = "param_" + int_fmt_str.format(param_ind)
1525 param_map[param] = param_str
1526 param_dct[param_str] = 0.0
1528 # Replace the parameters according to the map
1529 # Sort by string length (long to short) to prevent cases like x11
1530 # becoming f"{param_map["x1"]}1"
1531 for param in sorted(params, key=lambda s: -1.0 * len(s)):
1532 expression = expression.replace(param, param_map[param])
1534 # Partial linearity check
1535 for express_sec in expression.split("+"):
1536 in_sec = [param in express_sec for param in param_dct]
1537 n_params_in_sec = len(np.where(np.array(in_sec))[0])
1538 if n_params_in_sec > 1:
1539 raise ValueError(
1540 "FixParametricRelations expressions must be linear.")
1542 const_shift[expr_ind] = float(
1543 eval_expression(expression, param_dct))
1545 for param_ind in range(len(params)):
1546 param_str = "param_" + int_fmt_str.format(param_ind)
1547 if param_str not in expression:
1548 Jacobian[expr_ind, param_ind] = 0.0
1549 continue
1550 param_dct[param_str] = 1.0
1551 test_1 = float(eval_expression(expression, param_dct))
1552 test_1 -= const_shift[expr_ind]
1553 Jacobian[expr_ind, param_ind] = test_1
1555 param_dct[param_str] = 2.0
1556 test_2 = float(eval_expression(expression, param_dct))
1557 test_2 -= const_shift[expr_ind]
1558 if abs(test_2 / test_1 - 2.0) > eps:
1559 raise ValueError(
1560 "FixParametricRelations expressions must be linear.")
1561 param_dct[param_str] = 0.0
1563 args = [
1564 indices,
1565 Jacobian,
1566 const_shift,
1567 params,
1568 eps,
1569 use_cell,
1570 ]
1571 if cls is FixScaledParametricRelations:
1572 args = args[:-1]
1573 return cls(*args)
1575 @property
1576 def expressions(self):
1577 """Generate the expressions represented by the current self.Jacobian
1578 and self.const_shift objects"""
1579 expressions = []
1580 per = int(round(-1 * np.log10(self.eps)))
1581 fmt_str = "{:." + str(per + 1) + "g}"
1582 for index, shift_val in enumerate(self.const_shift):
1583 exp = ""
1584 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(
1585 shift_val) > self.eps:
1586 exp += fmt_str.format(shift_val)
1588 param_exp = ""
1589 for param_index, jacob_val in enumerate(self.Jacobian[index]):
1590 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)
1591 if abs_jacob_val < self.eps:
1592 continue
1594 param = self.params[param_index]
1595 if param_exp or exp:
1596 if jacob_val > -1.0 * self.eps:
1597 param_exp += " + "
1598 else:
1599 param_exp += " - "
1600 elif (not exp) and (not param_exp) and (
1601 jacob_val < -1.0 * self.eps):
1602 param_exp += "-"
1604 if np.abs(abs_jacob_val - 1.0) <= self.eps:
1605 param_exp += f"{param:s}"
1606 else:
1607 param_exp += (fmt_str +
1608 "*{:s}").format(abs_jacob_val, param)
1610 exp += param_exp
1612 expressions.append(exp)
1613 return np.array(expressions).reshape((-1, 3))
1615 def todict(self):
1616 """Create a dictionary representation of the constraint"""
1617 return {
1618 "name": type(self).__name__,
1619 "kwargs": {
1620 "indices": self.indices,
1621 "params": self.params,
1622 "Jacobian": self.Jacobian,
1623 "const_shift": self.const_shift,
1624 "eps": self.eps,
1625 "use_cell": self.use_cell,
1626 }
1627 }
1629 def __repr__(self):
1630 """The str representation of the constraint"""
1631 if len(self.indices) > 1:
1632 indices_str = "[{:d}, ..., {:d}]".format(
1633 self.indices[0], self.indices[-1])
1634 else:
1635 indices_str = f"[{self.indices[0]:d}]"
1637 if len(self.params) > 1:
1638 params_str = "[{:s}, ..., {:s}]".format(
1639 self.params[0], self.params[-1])
1640 elif len(self.params) == 1:
1641 params_str = f"[{self.params[0]:s}]"
1642 else:
1643 params_str = "[]"
1645 return '{:s}({:s}, {:s}, ..., {:e})'.format(
1646 type(self).__name__,
1647 indices_str,
1648 params_str,
1649 self.eps
1650 )
1653class FixScaledParametricRelations(FixParametricRelations):
1655 def __init__(
1656 self,
1657 indices,
1658 Jacobian,
1659 const_shift,
1660 params=None,
1661 eps=1e-12,
1662 ):
1663 """The fractional coordinate version of FixParametricRelations
1665 All arguments are the same, but since this is for fractional
1666 coordinates use_cell is false"""
1667 super().__init__(
1668 indices,
1669 Jacobian,
1670 const_shift,
1671 params,
1672 eps,
1673 False,
1674 )
1676 def adjust_contravariant(self, cell, vecs, B):
1677 """Adjust the values of a set of vectors that are contravariant
1678 with the unit transformation"""
1679 scaled = cell.scaled_positions(vecs).flatten()
1680 scaled = self.Jacobian_inv @ (scaled - B)
1681 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))
1683 return cell.cartesian_positions(scaled)
1685 def adjust_positions(self, atoms, positions):
1686 """Adjust positions of the atoms to match the constraints"""
1687 positions[self.indices] = self.adjust_contravariant(
1688 atoms.cell,
1689 positions[self.indices],
1690 self.const_shift,
1691 )
1692 positions[self.indices] = self.adjust_B(
1693 atoms.cell, positions[self.indices])
1695 def adjust_B(self, cell, positions):
1696 """Wraps the positions back to the unit cell and adjust B to
1697 keep track of this change"""
1698 fractional = cell.scaled_positions(positions)
1699 wrapped_fractional = (fractional % 1.0) % 1.0
1700 self.const_shift += np.round(wrapped_fractional - fractional).flatten()
1701 return cell.cartesian_positions(wrapped_fractional)
1703 def adjust_momenta(self, atoms, momenta):
1704 """Adjust momenta of the atoms to match the constraints"""
1705 momenta[self.indices] = self.adjust_contravariant(
1706 atoms.cell,
1707 momenta[self.indices],
1708 np.zeros(self.const_shift.shape),
1709 )
1711 def adjust_forces(self, atoms, forces):
1712 """Adjust forces of the atoms to match the constraints"""
1713 # Forces are coavarient to the coordinate transformation, use the
1714 # inverse transformations
1715 cart2frac_jacob = np.zeros(2 * (3 * len(atoms),))
1716 for i_atom in range(len(atoms)):
1717 cart2frac_jacob[3 * i_atom:3 * (i_atom + 1),
1718 3 * i_atom:3 * (i_atom + 1)] = atoms.cell.T
1720 jacobian = cart2frac_jacob @ self.Jacobian
1721 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T
1723 reduced_forces = jacobian.T @ forces.flatten()
1724 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)
1726 def todict(self):
1727 """Create a dictionary representation of the constraint"""
1728 dct = super().todict()
1729 del dct["kwargs"]["use_cell"]
1730 return dct
1733class FixCartesianParametricRelations(FixParametricRelations):
1735 def __init__(
1736 self,
1737 indices,
1738 Jacobian,
1739 const_shift,
1740 params=None,
1741 eps=1e-12,
1742 use_cell=False,
1743 ):
1744 """The Cartesian coordinate version of FixParametricRelations"""
1745 super().__init__(
1746 indices,
1747 Jacobian,
1748 const_shift,
1749 params,
1750 eps,
1751 use_cell,
1752 )
1754 def adjust_contravariant(self, vecs, B):
1755 """Adjust the values of a set of vectors that are contravariant with
1756 the unit transformation"""
1757 vecs = self.Jacobian_inv @ (vecs.flatten() - B)
1758 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))
1759 return vecs
1761 def adjust_positions(self, atoms, positions):
1762 """Adjust positions of the atoms to match the constraints"""
1763 if self.use_cell:
1764 return
1765 positions[self.indices] = self.adjust_contravariant(
1766 positions[self.indices],
1767 self.const_shift,
1768 )
1770 def adjust_momenta(self, atoms, momenta):
1771 """Adjust momenta of the atoms to match the constraints"""
1772 if self.use_cell:
1773 return
1774 momenta[self.indices] = self.adjust_contravariant(
1775 momenta[self.indices],
1776 np.zeros(self.const_shift.shape),
1777 )
1779 def adjust_forces(self, atoms, forces):
1780 """Adjust forces of the atoms to match the constraints"""
1781 if self.use_cell:
1782 return
1784 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()
1785 forces[self.indices] = (self.Jacobian_inv.T @
1786 forces_reduced).reshape(-1, 3)
1788 def adjust_cell(self, atoms, cell):
1789 """Adjust the cell of the atoms to match the constraints"""
1790 if not self.use_cell:
1791 return
1792 cell[self.indices] = self.adjust_contravariant(
1793 cell[self.indices],
1794 np.zeros(self.const_shift.shape),
1795 )
1797 def adjust_stress(self, atoms, stress):
1798 """Adjust the stress of the atoms to match the constraints"""
1799 if not self.use_cell:
1800 return
1802 stress_3x3 = voigt_6_to_full_3x3_stress(stress)
1803 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()
1804 stress_3x3[self.indices] = (
1805 self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)
1807 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)
1810class Hookean(FixConstraint):
1811 """Applies a Hookean restorative force between a pair of atoms, an atom
1812 and a point, or an atom and a plane."""
1814 def __init__(self, a1, a2, k, rt=None):
1815 """Forces two atoms to stay close together by applying no force if
1816 they are below a threshold length, rt, and applying a Hookean
1817 restorative force when the distance between them exceeds rt. Can
1818 also be used to tether an atom to a fixed point in space or to a
1819 distance above a plane.
1821 a1 : int
1822 Index of atom 1
1823 a2 : one of three options
1824 1) index of atom 2
1825 2) a fixed point in cartesian space to which to tether a1
1826 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.
1827 k : float
1828 Hooke's law (spring) constant to apply when distance
1829 exceeds threshold_length. Units of eV A^-2.
1830 rt : float
1831 The threshold length below which there is no force. The
1832 length is 1) between two atoms, 2) between atom and point.
1833 This argument is not supplied in case 3. Units of A.
1835 If a plane is specified, the Hooke's law force is applied if the atom
1836 is on the normal side of the plane. For instance, the plane with
1837 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z
1838 intercept of +7 and a normal vector pointing in the +z direction.
1839 If the atom has z > 7, then a downward force would be applied of
1840 k * (atom.z - 7). The same plane with the normal vector pointing in
1841 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).
1843 References:
1845 Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014)
1846 https://link.springer.com/article/10.1007%2Fs11244-013-0161-8
1847 """
1849 if isinstance(a2, int):
1850 self._type = 'two atoms'
1851 self.indices = [a1, a2]
1852 elif len(a2) == 3:
1853 self._type = 'point'
1854 self.index = a1
1855 self.origin = np.array(a2)
1856 elif len(a2) == 4:
1857 self._type = 'plane'
1858 self.index = a1
1859 self.plane = a2
1860 else:
1861 raise RuntimeError('Unknown type for a2')
1862 self.threshold = rt
1863 self.spring = k
1865 def get_removed_dof(self, atoms):
1866 return 0
1868 def todict(self):
1869 dct = {'name': 'Hookean'}
1870 dct['kwargs'] = {'rt': self.threshold,
1871 'k': self.spring}
1872 if self._type == 'two atoms':
1873 dct['kwargs']['a1'] = self.indices[0]
1874 dct['kwargs']['a2'] = self.indices[1]
1875 elif self._type == 'point':
1876 dct['kwargs']['a1'] = self.index
1877 dct['kwargs']['a2'] = self.origin
1878 elif self._type == 'plane':
1879 dct['kwargs']['a1'] = self.index
1880 dct['kwargs']['a2'] = self.plane
1881 else:
1882 raise NotImplementedError(f'Bad type: {self._type}')
1883 return dct
1885 def adjust_positions(self, atoms, newpositions):
1886 pass
1888 def adjust_momenta(self, atoms, momenta):
1889 pass
1891 def adjust_forces(self, atoms, forces):
1892 positions = atoms.positions
1893 if self._type == 'plane':
1894 A, B, C, D = self.plane
1895 x, y, z = positions[self.index]
1896 d = ((A * x + B * y + C * z + D) /
1897 np.sqrt(A**2 + B**2 + C**2))
1898 if d < 0:
1899 return
1900 magnitude = self.spring * d
1901 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))
1902 forces[self.index] += direction * magnitude
1903 return
1904 if self._type == 'two atoms':
1905 p1, p2 = positions[self.indices]
1906 elif self._type == 'point':
1907 p1 = positions[self.index]
1908 p2 = self.origin
1909 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1910 bondlength = np.linalg.norm(displace)
1911 if bondlength > self.threshold:
1912 magnitude = self.spring * (bondlength - self.threshold)
1913 direction = displace / np.linalg.norm(displace)
1914 if self._type == 'two atoms':
1915 forces[self.indices[0]] += direction * magnitude
1916 forces[self.indices[1]] -= direction * magnitude
1917 else:
1918 forces[self.index] += direction * magnitude
1920 def adjust_potential_energy(self, atoms):
1921 """Returns the difference to the potential energy due to an active
1922 constraint. (That is, the quantity returned is to be added to the
1923 potential energy.)"""
1924 positions = atoms.positions
1925 if self._type == 'plane':
1926 A, B, C, D = self.plane
1927 x, y, z = positions[self.index]
1928 d = ((A * x + B * y + C * z + D) /
1929 np.sqrt(A**2 + B**2 + C**2))
1930 if d > 0:
1931 return 0.5 * self.spring * d**2
1932 else:
1933 return 0.
1934 if self._type == 'two atoms':
1935 p1, p2 = positions[self.indices]
1936 elif self._type == 'point':
1937 p1 = positions[self.index]
1938 p2 = self.origin
1939 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1940 bondlength = np.linalg.norm(displace)
1941 if bondlength > self.threshold:
1942 return 0.5 * self.spring * (bondlength - self.threshold)**2
1943 else:
1944 return 0.
1946 def get_indices(self):
1947 if self._type == 'two atoms':
1948 return self.indices
1949 elif self._type == 'point':
1950 return self.index
1951 elif self._type == 'plane':
1952 return self.index
1954 def index_shuffle(self, atoms, ind):
1955 # See docstring of superclass
1956 if self._type == 'two atoms':
1957 newa = [-1, -1] # Signal error
1958 for new, old in slice2enlist(ind, len(atoms)):
1959 for i, a in enumerate(self.indices):
1960 if old == a:
1961 newa[i] = new
1962 if newa[0] == -1 or newa[1] == -1:
1963 raise IndexError('Constraint not part of slice')
1964 self.indices = newa
1965 elif (self._type == 'point') or (self._type == 'plane'):
1966 newa = -1 # Signal error
1967 for new, old in slice2enlist(ind, len(atoms)):
1968 if old == self.index:
1969 newa = new
1970 break
1971 if newa == -1:
1972 raise IndexError('Constraint not part of slice')
1973 self.index = newa
1975 def __repr__(self):
1976 if self._type == 'two atoms':
1977 return 'Hookean(%d, %d)' % tuple(self.indices)
1978 elif self._type == 'point':
1979 return 'Hookean(%d) to cartesian' % self.index
1980 else:
1981 return 'Hookean(%d) to plane' % self.index
1984class ExternalForce(FixConstraint):
1985 """Constraint object for pulling two atoms apart by an external force.
1987 You can combine this constraint for example with FixBondLength but make
1988 sure that *ExternalForce* comes first in the list if there are overlaps
1989 between atom1-2 and atom3-4:
1991 >>> from ase.build import bulk
1993 >>> atoms = bulk('Cu', 'fcc', a=3.6)
1994 >>> atom1, atom2, atom3, atom4 = atoms[:4]
1995 >>> fext = 1.0
1996 >>> con1 = ExternalForce(atom1, atom2, f_ext)
1997 >>> con2 = FixBondLength(atom3, atom4)
1998 >>> atoms.set_constraint([con1, con2])
2000 see ase/test/external_force.py"""
2002 def __init__(self, a1, a2, f_ext):
2003 self.indices = [a1, a2]
2004 self.external_force = f_ext
2006 def get_removed_dof(self, atoms):
2007 return 0
2009 def adjust_positions(self, atoms, new):
2010 pass
2012 def adjust_forces(self, atoms, forces):
2013 dist = np.subtract.reduce(atoms.positions[self.indices])
2014 force = self.external_force * dist / np.linalg.norm(dist)
2015 forces[self.indices] += (force, -force)
2017 def adjust_potential_energy(self, atoms):
2018 dist = np.subtract.reduce(atoms.positions[self.indices])
2019 return -np.linalg.norm(dist) * self.external_force
2021 def index_shuffle(self, atoms, ind):
2022 """Shuffle the indices of the two atoms in this constraint"""
2023 newa = [-1, -1] # Signal error
2024 for new, old in slice2enlist(ind, len(atoms)):
2025 for i, a in enumerate(self.indices):
2026 if old == a:
2027 newa[i] = new
2028 if newa[0] == -1 or newa[1] == -1:
2029 raise IndexError('Constraint not part of slice')
2030 self.indices = newa
2032 def __repr__(self):
2033 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],
2034 self.indices[1],
2035 self.external_force)
2037 def todict(self):
2038 return {'name': 'ExternalForce',
2039 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2040 'f_ext': self.external_force}}
2043class MirrorForce(FixConstraint):
2044 """Constraint object for mirroring the force between two atoms.
2046 This class is designed to find a transition state with the help of a
2047 single optimization. It can be used if the transition state belongs to a
2048 bond breaking reaction. First the given bond length will be fixed until
2049 all other degrees of freedom are optimized, then the forces of the two
2050 atoms will be mirrored to find the transition state. The mirror plane is
2051 perpendicular to the connecting line of the atoms. Transition states in
2052 dependence of the force can be obtained by stretching the molecule and
2053 fixing its total length with *FixBondLength* or by using *ExternalForce*
2054 during the optimization with *MirrorForce*.
2056 Parameters
2057 ----------
2058 a1: int
2059 First atom index.
2060 a2: int
2061 Second atom index.
2062 max_dist: float
2063 Upper limit of the bond length interval where the transition state
2064 can be found.
2065 min_dist: float
2066 Lower limit of the bond length interval where the transition state
2067 can be found.
2068 fmax: float
2069 Maximum force used for the optimization.
2071 Notes
2072 -----
2073 You can combine this constraint for example with FixBondLength but make
2074 sure that *MirrorForce* comes first in the list if there are overlaps
2075 between atom1-2 and atom3-4:
2077 >>> from ase.build import bulk
2079 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2080 >>> atom1, atom2, atom3, atom4 = atoms[:4]
2081 >>> con1 = MirrorForce(atom1, atom2)
2082 >>> con2 = FixBondLength(atom3, atom4)
2083 >>> atoms.set_constraint([con1, con2])
2085 """
2087 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):
2088 self.indices = [a1, a2]
2089 self.min_dist = min_dist
2090 self.max_dist = max_dist
2091 self.fmax = fmax
2093 def adjust_positions(self, atoms, new):
2094 pass
2096 def adjust_forces(self, atoms, forces):
2097 dist = np.subtract.reduce(atoms.positions[self.indices])
2098 d = np.linalg.norm(dist)
2099 if (d < self.min_dist) or (d > self.max_dist):
2100 # Stop structure optimization
2101 forces[:] *= 0
2102 return
2103 dist /= d
2104 df = np.subtract.reduce(forces[self.indices])
2105 f = df.dot(dist)
2106 con_saved = atoms.constraints
2107 try:
2108 con = [con for con in con_saved
2109 if not isinstance(con, MirrorForce)]
2110 atoms.set_constraint(con)
2111 forces_copy = atoms.get_forces()
2112 finally:
2113 atoms.set_constraint(con_saved)
2114 df1 = -1 / 2. * f * dist
2115 forces_copy[self.indices] += (df1, -df1)
2116 # Check if forces would be converged if the bond with mirrored forces
2117 # would also be fixed
2118 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2119 factor = 1.
2120 else:
2121 factor = 0.
2122 df1 = -(1 + factor) / 2. * f * dist
2123 forces[self.indices] += (df1, -df1)
2125 def index_shuffle(self, atoms, ind):
2126 """Shuffle the indices of the two atoms in this constraint
2128 """
2129 newa = [-1, -1] # Signal error
2130 for new, old in slice2enlist(ind, len(atoms)):
2131 for i, a in enumerate(self.indices):
2132 if old == a:
2133 newa[i] = new
2134 if newa[0] == -1 or newa[1] == -1:
2135 raise IndexError('Constraint not part of slice')
2136 self.indices = newa
2138 def __repr__(self):
2139 return 'MirrorForce(%d, %d, %f, %f, %f)' % (
2140 self.indices[0], self.indices[1], self.max_dist, self.min_dist,
2141 self.fmax)
2143 def todict(self):
2144 return {'name': 'MirrorForce',
2145 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2146 'max_dist': self.max_dist,
2147 'min_dist': self.min_dist, 'fmax': self.fmax}}
2150class MirrorTorque(FixConstraint):
2151 """Constraint object for mirroring the torque acting on a dihedral
2152 angle defined by four atoms.
2154 This class is designed to find a transition state with the help of a
2155 single optimization. It can be used if the transition state belongs to a
2156 cis-trans-isomerization with a change of dihedral angle. First the given
2157 dihedral angle will be fixed until all other degrees of freedom are
2158 optimized, then the torque acting on the dihedral angle will be mirrored
2159 to find the transition state. Transition states in
2160 dependence of the force can be obtained by stretching the molecule and
2161 fixing its total length with *FixBondLength* or by using *ExternalForce*
2162 during the optimization with *MirrorTorque*.
2164 This constraint can be used to find
2165 transition states of cis-trans-isomerization.
2167 a1 a4
2168 | |
2169 a2 __ a3
2171 Parameters
2172 ----------
2173 a1: int
2174 First atom index.
2175 a2: int
2176 Second atom index.
2177 a3: int
2178 Third atom index.
2179 a4: int
2180 Fourth atom index.
2181 max_angle: float
2182 Upper limit of the dihedral angle interval where the transition state
2183 can be found.
2184 min_angle: float
2185 Lower limit of the dihedral angle interval where the transition state
2186 can be found.
2187 fmax: float
2188 Maximum force used for the optimization.
2190 Notes
2191 -----
2192 You can combine this constraint for example with FixBondLength but make
2193 sure that *MirrorTorque* comes first in the list if there are overlaps
2194 between atom1-4 and atom5-6:
2196 >>> from ase.build import bulk
2198 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2199 >>> atom1, atom2, atom3, atom4, atom5, atom6 = atoms[:6]
2200 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)
2201 >>> con2 = FixBondLength(atom5, atom6)
2202 >>> atoms.set_constraint([con1, con2])
2204 """
2206 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,
2207 fmax=0.1):
2208 self.indices = [a1, a2, a3, a4]
2209 self.min_angle = min_angle
2210 self.max_angle = max_angle
2211 self.fmax = fmax
2213 def adjust_positions(self, atoms, new):
2214 pass
2216 def adjust_forces(self, atoms, forces):
2217 angle = atoms.get_dihedral(self.indices[0], self.indices[1],
2218 self.indices[2], self.indices[3])
2219 angle *= np.pi / 180.
2220 if (angle < self.min_angle) or (angle > self.max_angle):
2221 # Stop structure optimization
2222 forces[:] *= 0
2223 return
2224 p = atoms.positions[self.indices]
2225 f = forces[self.indices]
2227 f0 = (f[1] + f[2]) / 2.
2228 ff = f - f0
2229 p0 = (p[2] + p[1]) / 2.
2230 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)
2231 fff = ff - np.cross(m0, p - p0)
2232 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \
2233 (p[1] - p0).dot(p[1] - p0)
2234 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \
2235 (p[2] - p0).dot(p[2] - p0)
2236 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \
2237 np.linalg.norm(p[1] - p0)
2238 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \
2239 np.linalg.norm(p[2] - p0)
2240 omegap = omegap1 + omegap2
2241 con_saved = atoms.constraints
2242 try:
2243 con = [con for con in con_saved
2244 if not isinstance(con, MirrorTorque)]
2245 atoms.set_constraint(con)
2246 forces_copy = atoms.get_forces()
2247 finally:
2248 atoms.set_constraint(con_saved)
2249 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \
2250 np.linalg.norm(p[1] - p0)
2251 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \
2252 np.linalg.norm(p[2] - p0)
2253 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2254 # Check if forces would be converged if the dihedral angle with
2255 # mirrored torque would also be fixed
2256 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2257 factor = 1.
2258 else:
2259 factor = 0.
2260 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \
2261 np.linalg.norm(p[1] - p0)
2262 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \
2263 np.linalg.norm(p[2] - p0)
2264 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2266 def index_shuffle(self, atoms, ind):
2267 # See docstring of superclass
2268 indices = []
2269 for new, old in slice2enlist(ind, len(atoms)):
2270 if old in self.indices:
2271 indices.append(new)
2272 if len(indices) == 0:
2273 raise IndexError('All indices in MirrorTorque not part of slice')
2274 self.indices = np.asarray(indices, int)
2276 def __repr__(self):
2277 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (
2278 self.indices[0], self.indices[1], self.indices[2],
2279 self.indices[3], self.max_angle, self.min_angle, self.fmax)
2281 def todict(self):
2282 return {'name': 'MirrorTorque',
2283 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2284 'a3': self.indices[2], 'a4': self.indices[3],
2285 'max_angle': self.max_angle,
2286 'min_angle': self.min_angle, 'fmax': self.fmax}}
2289class FixSymmetry(FixConstraint):
2290 """
2291 Constraint to preserve spacegroup symmetry during optimisation.
2293 Requires spglib package to be available.
2294 """
2296 def __init__(self, atoms, symprec=0.01, adjust_positions=True,
2297 adjust_cell=True, verbose=False):
2298 self.atoms = atoms.copy()
2299 self.symprec = symprec
2300 self.verbose = verbose
2301 refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry
2302 sym = prep_symmetry(atoms, symprec, self.verbose)
2303 self.rotations, self.translations, self.symm_map = sym
2304 self.do_adjust_positions = adjust_positions
2305 self.do_adjust_cell = adjust_cell
2307 def adjust_cell(self, atoms, cell):
2308 if not self.do_adjust_cell:
2309 return
2310 # stress should definitely be symmetrized as a rank 2 tensor
2311 # UnitCellFilter uses deformation gradient as cell DOF with steps
2312 # dF = stress.F^-T quantity that should be symmetrized is therefore dF .
2313 # F^T assume prev F = I, so just symmetrize dF
2314 cur_cell = atoms.get_cell()
2315 cur_cell_inv = atoms.cell.reciprocal().T
2317 # F defined such that cell = cur_cell . F^T
2318 # assume prev F = I, so dF = F - I
2319 delta_deform_grad = np.dot(cur_cell_inv, cell).T - np.eye(3)
2321 # symmetrization doesn't work properly with large steps, since
2322 # it depends on current cell, and cell is being changed by deformation
2323 # gradient
2324 max_delta_deform_grad = np.max(np.abs(delta_deform_grad))
2325 if max_delta_deform_grad > 0.25:
2326 raise RuntimeError('FixSymmetry adjust_cell does not work properly'
2327 ' with large deformation gradient step {} > 0.25'
2328 .format(max_delta_deform_grad))
2329 elif max_delta_deform_grad > 0.15:
2330 warn('FixSymmetry adjust_cell may be ill behaved with large '
2331 'deformation gradient step {}'.format(max_delta_deform_grad))
2333 symmetrized_delta_deform_grad = symmetrize_rank2(cur_cell, cur_cell_inv,
2334 delta_deform_grad,
2335 self.rotations)
2336 cell[:] = np.dot(cur_cell,
2337 (symmetrized_delta_deform_grad + np.eye(3)).T)
2339 def adjust_positions(self, atoms, new):
2340 if not self.do_adjust_positions:
2341 return
2342 # symmetrize changes in position as rank 1 tensors
2343 step = new - atoms.positions
2344 symmetrized_step = symmetrize_rank1(atoms.get_cell(),
2345 atoms.cell.reciprocal().T, step,
2346 self.rotations, self.translations,
2347 self.symm_map)
2348 new[:] = atoms.positions + symmetrized_step
2350 def adjust_forces(self, atoms, forces):
2351 # symmetrize forces as rank 1 tensors
2352 # print('adjusting forces')
2353 forces[:] = symmetrize_rank1(atoms.get_cell(),
2354 atoms.cell.reciprocal().T, forces,
2355 self.rotations, self.translations,
2356 self.symm_map)
2358 def adjust_stress(self, atoms, stress):
2359 # symmetrize stress as rank 2 tensor
2360 raw_stress = voigt_6_to_full_3x3_stress(stress)
2361 symmetrized_stress = symmetrize_rank2(atoms.get_cell(),
2362 atoms.cell.reciprocal().T,
2363 raw_stress, self.rotations)
2364 stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress)
2366 def index_shuffle(self, atoms, ind):
2367 if len(atoms) != len(ind) or len(set(ind)) != len(ind):
2368 raise RuntimeError("FixSymmetry can only accomodate atom"
2369 " permutions, and len(Atoms) == {} "
2370 "!= len(ind) == {} or ind has duplicates"
2371 .format(len(atoms), len(ind)))
2373 ind_reversed = np.zeros((len(ind)), dtype=int)
2374 ind_reversed[ind] = range(len(ind))
2375 new_symm_map = []
2376 for sm in self.symm_map:
2377 new_sm = np.array([-1] * len(atoms))
2378 for at_i in range(len(ind)):
2379 new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]]
2380 new_symm_map.append(new_sm)
2382 self.symm_map = new_symm_map
2384 def todict(self):
2385 return {
2386 'name': 'FixSymmetry',
2387 'kwargs': {
2388 'atoms': self.atoms,
2389 'symprec': self.symprec,
2390 'adjust_positions': self.do_adjust_positions,
2391 'adjust_cell': self.do_adjust_cell,
2392 'verbose': self.verbose,
2393 },
2394 }
2397class Filter(FilterOld):
2398 @deprecated('Import Filter from ase.filters')
2399 def __init__(self, *args, **kwargs):
2400 """
2401 .. deprecated:: 3.23.0
2402 Import ``Filter`` from :mod:`ase.filters`
2403 """
2404 super().__init__(*args, **kwargs)
2407class StrainFilter(StrainFilterOld):
2408 @deprecated('Import StrainFilter from ase.filters')
2409 def __init__(self, *args, **kwargs):
2410 """
2411 .. deprecated:: 3.23.0
2412 Import ``StrainFilter`` from :mod:`ase.filters`
2413 """
2414 super().__init__(*args, **kwargs)
2417class UnitCellFilter(UnitCellFilterOld):
2418 @deprecated('Import UnitCellFilter from ase.filters')
2419 def __init__(self, *args, **kwargs):
2420 """
2421 .. deprecated:: 3.23.0
2422 Import ``UnitCellFilter`` from :mod:`ase.filters`
2423 """
2424 super().__init__(*args, **kwargs)
2427class ExpCellFilter(ExpCellFilterOld):
2428 @deprecated('Import ExpCellFilter from ase.filters')
2429 def __init__(self, *args, **kwargs):
2430 """
2431 .. deprecated:: 3.23.0
2432 Import ``ExpCellFilter`` from :mod:`ase.filters`
2433 or use :class:`~ase.filters.FrechetCellFilter` for better
2434 convergence w.r.t. cell variables
2435 """
2436 super().__init__(*args, **kwargs)