Coverage for /builds/kinetik161/ase/ase/constraints.py: 86.76%
1156 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Constraints"""
2from typing import Sequence
3from warnings import warn
5import numpy as np
7from ase.filters import ExpCellFilter as ExpCellFilterOld
8from ase.filters import Filter as FilterOld
9from ase.filters import StrainFilter as StrainFilterOld
10from ase.filters import UnitCellFilter as UnitCellFilterOld
11from ase.geometry import (conditional_find_mic, find_mic, get_angles,
12 get_angles_derivatives, get_dihedrals,
13 get_dihedrals_derivatives, get_distances_derivatives,
14 wrap_positions)
15from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
16from ase.utils import deprecated
17from ase.utils.parsemath import eval_expression
19__all__ = [
20 'FixCartesian', 'FixBondLength', 'FixedMode',
21 'FixAtoms', 'FixScaled', 'FixCom', 'FixedPlane',
22 'FixConstraint', 'FixedLine', 'FixBondLengths', 'FixLinearTriatomic',
23 'FixInternals', 'Hookean', 'ExternalForce', 'MirrorForce', 'MirrorTorque',
24 "FixScaledParametricRelations", "FixCartesianParametricRelations"]
27def dict2constraint(dct):
28 if dct['name'] not in __all__:
29 raise ValueError
30 return globals()[dct['name']](**dct['kwargs'])
33def slice2enlist(s, n):
34 """Convert a slice object into a list of (new, old) tuples."""
35 if isinstance(s, slice):
36 return enumerate(range(*s.indices(n)))
37 return enumerate(s)
40def constrained_indices(atoms, only_include=None):
41 """Returns a list of indices for the atoms that are constrained
42 by a constraint that is applied. By setting only_include to a
43 specific type of constraint you can make it only look for that
44 given constraint.
45 """
46 indices = []
47 for constraint in atoms.constraints:
48 if only_include is not None:
49 if not isinstance(constraint, only_include):
50 continue
51 indices.extend(np.array(constraint.get_indices()))
52 return np.array(np.unique(indices))
55class FixConstraint:
56 """Base class for classes that fix one or more atoms in some way."""
58 def index_shuffle(self, atoms, ind):
59 """Change the indices.
61 When the ordering of the atoms in the Atoms object changes,
62 this method can be called to shuffle the indices of the
63 constraints.
65 ind -- List or tuple of indices.
67 """
68 raise NotImplementedError
70 def repeat(self, m, n):
71 """ basic method to multiply by m, needs to know the length
72 of the underlying atoms object for the assignment of
73 multiplied constraints to work.
74 """
75 msg = ("Repeat is not compatible with your atoms' constraints."
76 ' Use atoms.set_constraint() before calling repeat to '
77 'remove your constraints.')
78 raise NotImplementedError(msg)
80 def adjust_momenta(self, atoms, momenta):
81 """Adjusts momenta in identical manner to forces."""
82 self.adjust_forces(atoms, momenta)
84 def copy(self):
85 return dict2constraint(self.todict().copy())
88class IndexedConstraint(FixConstraint):
89 def __init__(self, indices=None, mask=None):
90 """Constrain chosen atoms.
92 Parameters
93 ----------
94 indices : sequence of int
95 Indices for those atoms that should be constrained.
96 mask : sequence of bool
97 One boolean per atom indicating if the atom should be
98 constrained or not.
99 """
101 if mask is not None:
102 if indices is not None:
103 raise ValueError('Use only one of "indices" and "mask".')
104 indices = mask
105 indices = np.atleast_1d(indices)
106 if np.ndim(indices) > 1:
107 raise ValueError('indices has wrong amount of dimensions. '
108 f'Got {np.ndim(indices)}, expected ndim <= 1')
110 if indices.dtype == bool:
111 indices = np.arange(len(indices))[indices]
112 elif len(indices) == 0:
113 indices = np.empty(0, dtype=int)
114 elif not np.issubdtype(indices.dtype, np.integer):
115 raise ValueError('Indices must be integers or boolean mask, '
116 f'not dtype={indices.dtype}')
118 if len(set(indices)) < len(indices):
119 raise ValueError(
120 'The indices array contains duplicates. '
121 'Perhaps you want to specify a mask instead, but '
122 'forgot the mask= keyword.')
124 self.index = indices
126 def index_shuffle(self, atoms, ind):
127 # See docstring of superclass
128 index = []
130 # Resolve negative indices:
131 actual_indices = set(np.arange(len(atoms))[self.index])
133 for new, old in slice2enlist(ind, len(atoms)):
134 if old in actual_indices:
135 index.append(new)
136 if len(index) == 0:
137 raise IndexError('All indices in FixAtoms not part of slice')
138 self.index = np.asarray(index, int)
139 # XXX make immutable
141 def get_indices(self):
142 return self.index.copy()
144 def repeat(self, m, n):
145 i0 = 0
146 natoms = 0
147 if isinstance(m, int):
148 m = (m, m, m)
149 index_new = []
150 for m2 in range(m[2]):
151 for m1 in range(m[1]):
152 for m0 in range(m[0]):
153 i1 = i0 + n
154 index_new += [i + natoms for i in self.index]
155 i0 = i1
156 natoms += n
157 self.index = np.asarray(index_new, int)
158 # XXX make immutable
159 return self
161 def delete_atoms(self, indices, natoms):
162 """Removes atoms from the index array, if present.
164 Required for removing atoms with existing constraint.
165 """
167 i = np.zeros(natoms, int) - 1
168 new = np.delete(np.arange(natoms), indices)
169 i[new] = np.arange(len(new))
170 index = i[self.index]
171 self.index = index[index >= 0]
172 # XXX make immutable
173 if len(self.index) == 0:
174 return None
175 return self
178class FixAtoms(IndexedConstraint):
179 """Fix chosen atoms.
181 Examples
182 --------
183 Fix all Copper atoms:
185 >>> from ase.build import bulk
187 >>> atoms = bulk('Cu', 'fcc', a=3.6)
188 >>> mask = (atoms.symbols == 'Cu')
189 >>> c = FixAtoms(mask=mask)
190 >>> atoms.set_constraint(c)
192 Fix all atoms with z-coordinate less than 1.0 Angstrom:
194 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)
195 >>> atoms.set_constraint(c)
196 """
198 def get_removed_dof(self, atoms):
199 return 3 * len(self.index)
201 def adjust_positions(self, atoms, new):
202 new[self.index] = atoms.positions[self.index]
204 def adjust_forces(self, atoms, forces):
205 forces[self.index] = 0.0
207 def __repr__(self):
208 clsname = type(self).__name__
209 indices = ints2string(self.index)
210 return f'{clsname}(indices={indices})'
212 def todict(self):
213 return {'name': 'FixAtoms',
214 'kwargs': {'indices': self.index.tolist()}}
217class FixCom(FixConstraint):
218 """Constraint class for fixing the center of mass."""
220 def get_removed_dof(self, atoms):
221 return 3
223 def adjust_positions(self, atoms, new):
224 masses = atoms.get_masses()
225 old_cm = atoms.get_center_of_mass()
226 new_cm = np.dot(masses, new) / masses.sum()
227 diff = old_cm - new_cm
228 new += diff
230 def adjust_momenta(self, atoms, momenta):
231 """Adjust momenta so that the center-of-mass velocity is zero."""
232 masses = atoms.get_masses()
233 velocity_com = momenta.sum(axis=0) / masses.sum()
234 momenta -= masses[:, None] * velocity_com
236 def adjust_forces(self, atoms, forces):
237 # Eqs. (3) and (7) in https://doi.org/10.1021/jp9722824
238 masses = atoms.get_masses()
239 forces -= masses[:, None] * (masses @ forces) / sum(masses**2)
241 def todict(self):
242 return {'name': 'FixCom',
243 'kwargs': {}}
246def ints2string(x, threshold=None):
247 """Convert ndarray of ints to string."""
248 if threshold is None or len(x) <= threshold:
249 return str(x.tolist())
250 return str(x[:threshold].tolist())[:-1] + ', ...]'
253class FixBondLengths(FixConstraint):
254 maxiter = 500
256 def __init__(self, pairs, tolerance=1e-13,
257 bondlengths=None, iterations=None):
258 """iterations:
259 Ignored"""
260 self.pairs = np.asarray(pairs)
261 self.tolerance = tolerance
262 self.bondlengths = bondlengths
264 def get_removed_dof(self, atoms):
265 return len(self.pairs)
267 def adjust_positions(self, atoms, new):
268 old = atoms.positions
269 masses = atoms.get_masses()
271 if self.bondlengths is None:
272 self.bondlengths = self.initialize_bond_lengths(atoms)
274 for i in range(self.maxiter):
275 converged = True
276 for j, ab in enumerate(self.pairs):
277 a = ab[0]
278 b = ab[1]
279 cd = self.bondlengths[j]
280 r0 = old[a] - old[b]
281 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
282 d1 = new[a] - new[b] - r0 + d0
283 m = 1 / (1 / masses[a] + 1 / masses[b])
284 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)
285 if abs(x) > self.tolerance:
286 new[a] += x * m / masses[a] * d0
287 new[b] -= x * m / masses[b] * d0
288 converged = False
289 if converged:
290 break
291 else:
292 raise RuntimeError('Did not converge')
294 def adjust_momenta(self, atoms, p):
295 old = atoms.positions
296 masses = atoms.get_masses()
298 if self.bondlengths is None:
299 self.bondlengths = self.initialize_bond_lengths(atoms)
301 for i in range(self.maxiter):
302 converged = True
303 for j, ab in enumerate(self.pairs):
304 a = ab[0]
305 b = ab[1]
306 cd = self.bondlengths[j]
307 d = old[a] - old[b]
308 d, _ = find_mic(d, atoms.cell, atoms.pbc)
309 dv = p[a] / masses[a] - p[b] / masses[b]
310 m = 1 / (1 / masses[a] + 1 / masses[b])
311 x = -np.dot(dv, d) / cd**2
312 if abs(x) > self.tolerance:
313 p[a] += x * m * d
314 p[b] -= x * m * d
315 converged = False
316 if converged:
317 break
318 else:
319 raise RuntimeError('Did not converge')
321 def adjust_forces(self, atoms, forces):
322 self.constraint_forces = -forces
323 self.adjust_momenta(atoms, forces)
324 self.constraint_forces += forces
326 def initialize_bond_lengths(self, atoms):
327 bondlengths = np.zeros(len(self.pairs))
329 for i, ab in enumerate(self.pairs):
330 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)
332 return bondlengths
334 def get_indices(self):
335 return np.unique(self.pairs.ravel())
337 def todict(self):
338 return {'name': 'FixBondLengths',
339 'kwargs': {'pairs': self.pairs.tolist(),
340 'tolerance': self.tolerance}}
342 def index_shuffle(self, atoms, ind):
343 """Shuffle the indices of the two atoms in this constraint"""
344 map = np.zeros(len(atoms), int)
345 map[ind] = 1
346 n = map.sum()
347 map[:] = -1
348 map[ind] = range(n)
349 pairs = map[self.pairs]
350 self.pairs = pairs[(pairs != -1).all(1)]
351 if len(self.pairs) == 0:
352 raise IndexError('Constraint not part of slice')
355def FixBondLength(a1, a2):
356 """Fix distance between atoms with indices a1 and a2."""
357 return FixBondLengths([(a1, a2)])
360class FixLinearTriatomic(FixConstraint):
361 """Holonomic constraints for rigid linear triatomic molecules."""
363 def __init__(self, triples):
364 """Apply RATTLE-type bond constraints between outer atoms n and m
365 and linear vectorial constraints to the position of central
366 atoms o to fix the geometry of linear triatomic molecules of the
367 type:
369 n--o--m
371 Parameters:
373 triples: list
374 Indices of the atoms forming the linear molecules to constrain
375 as triples. Sequence should be (n, o, m) or (m, o, n).
377 When using these constraints in molecular dynamics or structure
378 optimizations, atomic forces need to be redistributed within a
379 triple. The function redistribute_forces_optimization implements
380 the redistribution of forces for structure optimization, while
381 the function redistribute_forces_md implements the redistribution
382 for molecular dynamics.
384 References:
386 Ciccotti et al. Molecular Physics 47 (1982)
387 :doi:`10.1080/00268978200100942`
388 """
389 self.triples = np.asarray(triples)
390 if self.triples.shape[1] != 3:
391 raise ValueError('"triples" has wrong size')
392 self.bondlengths = None
394 def get_removed_dof(self, atoms):
395 return 4 * len(self.triples)
397 @property
398 def n_ind(self):
399 return self.triples[:, 0]
401 @property
402 def m_ind(self):
403 return self.triples[:, 2]
405 @property
406 def o_ind(self):
407 return self.triples[:, 1]
409 def initialize(self, atoms):
410 masses = atoms.get_masses()
411 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)
413 self.bondlengths = self.initialize_bond_lengths(atoms)
414 self.bondlengths_nm = self.bondlengths.sum(axis=1)
416 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]
417 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +
418 C1[:, 1] ** 2 * self.mass_n * self.mass_o +
419 self.mass_n * self.mass_m)
420 C2 = C1 / C2[:, None]
421 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]
422 C3 = C2 * self.mass_o[:, None] * C3[:, None]
423 C3[:, 1] *= -1
424 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T
425 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)
426 C4 = C1 / C4[:, None]
428 self.C1 = C1
429 self.C2 = C2
430 self.C3 = C3
431 self.C4 = C4
433 def adjust_positions(self, atoms, new):
434 old = atoms.positions
435 new_n, new_m, new_o = self.get_slices(new)
437 if self.bondlengths is None:
438 self.initialize(atoms)
440 r0 = old[self.n_ind] - old[self.m_ind]
441 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
442 d1 = new_n - new_m - r0 + d0
443 a = np.einsum('ij,ij->i', d0, d0)
444 b = np.einsum('ij,ij->i', d1, d0)
445 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2
446 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))
447 g = g[:, None] * self.C3
448 new_n -= g[:, 0, None] * d0
449 new_m += g[:, 1, None] * d0
450 if np.allclose(d0, r0):
451 new_o = (self.C1[:, 0, None] * new_n
452 + self.C1[:, 1, None] * new_m)
453 else:
454 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)
455 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)
456 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2
457 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)
459 self.set_slices(new_n, new_m, new_o, new)
461 def adjust_momenta(self, atoms, p):
462 old = atoms.positions
463 p_n, p_m, p_o = self.get_slices(p)
465 if self.bondlengths is None:
466 self.initialize(atoms)
468 mass_nn = self.mass_n[:, None]
469 mass_mm = self.mass_m[:, None]
470 mass_oo = self.mass_o[:, None]
472 d = old[self.n_ind] - old[self.m_ind]
473 d, _ = find_mic(d, atoms.cell, atoms.pbc)
474 dv = p_n / mass_nn - p_m / mass_mm
475 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2
476 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]
477 p_n -= k[:, 0, None] * mass_nn * d
478 p_m += k[:, 1, None] * mass_mm * d
479 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +
480 self.C1[:, 1, None] * p_m / mass_mm))
482 self.set_slices(p_n, p_m, p_o, p)
484 def adjust_forces(self, atoms, forces):
486 if self.bondlengths is None:
487 self.initialize(atoms)
489 A = self.C4 * np.diff(self.C1)
490 A[:, 0] *= -1
491 A -= 1
492 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]
493 A /= (A.sum(axis=1))[:, None]
495 self.constraint_forces = -forces
496 old = atoms.positions
498 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)
500 d = old[self.n_ind] - old[self.m_ind]
501 d, _ = find_mic(d, atoms.cell, atoms.pbc)
502 df = fr_n - fr_m
503 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2
504 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]
505 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]
506 forces[self.o_ind] = fr_o + k[:, None] * d * B
508 self.constraint_forces += forces
510 def redistribute_forces_optimization(self, forces):
511 """Redistribute forces within a triple when performing structure
512 optimizations.
514 The redistributed forces needs to be further adjusted using the
515 appropriate Lagrange multipliers as implemented in adjust_forces."""
516 forces_n, forces_m, forces_o = self.get_slices(forces)
517 C1_1 = self.C1[:, 0, None]
518 C1_2 = self.C1[:, 1, None]
519 C4_1 = self.C4[:, 0, None]
520 C4_2 = self.C4[:, 1, None]
522 fr_n = ((1 - C4_1 * C1_1) * forces_n -
523 C4_1 * (C1_2 * forces_m - forces_o))
524 fr_m = ((1 - C4_2 * C1_2) * forces_m -
525 C4_2 * (C1_1 * forces_n - forces_o))
526 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +
527 C4_1 * forces_n + C4_2 * forces_m)
529 return fr_n, fr_m, fr_o
531 def redistribute_forces_md(self, atoms, forces, rand=False):
532 """Redistribute forces within a triple when performing molecular
533 dynamics.
535 When rand=True, use the equations for random force terms, as
536 used e.g. by Langevin dynamics, otherwise apply the standard
537 equations for deterministic forces (see Ciccotti et al. Molecular
538 Physics 47 (1982))."""
539 if self.bondlengths is None:
540 self.initialize(atoms)
541 forces_n, forces_m, forces_o = self.get_slices(forces)
542 C1_1 = self.C1[:, 0, None]
543 C1_2 = self.C1[:, 1, None]
544 C2_1 = self.C2[:, 0, None]
545 C2_2 = self.C2[:, 1, None]
546 mass_nn = self.mass_n[:, None]
547 mass_mm = self.mass_m[:, None]
548 mass_oo = self.mass_o[:, None]
549 if rand:
550 mr1 = (mass_mm / mass_nn) ** 0.5
551 mr2 = (mass_oo / mass_nn) ** 0.5
552 mr3 = (mass_nn / mass_mm) ** 0.5
553 mr4 = (mass_oo / mass_mm) ** 0.5
554 else:
555 mr1 = 1.0
556 mr2 = 1.0
557 mr3 = 1.0
558 mr4 = 1.0
560 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -
561 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -
562 mr2 * mass_mm * mass_nn * forces_o))
564 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -
565 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -
566 mr4 * mass_mm * mass_nn * forces_o))
568 self.set_slices(fr_n, fr_m, 0.0, forces)
570 def get_slices(self, a):
571 a_n = a[self.n_ind]
572 a_m = a[self.m_ind]
573 a_o = a[self.o_ind]
575 return a_n, a_m, a_o
577 def set_slices(self, a_n, a_m, a_o, a):
578 a[self.n_ind] = a_n
579 a[self.m_ind] = a_m
580 a[self.o_ind] = a_o
582 def initialize_bond_lengths(self, atoms):
583 bondlengths = np.zeros((len(self.triples), 2))
585 for i in range(len(self.triples)):
586 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],
587 self.o_ind[i], mic=True)
588 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],
589 self.m_ind[i], mic=True)
591 return bondlengths
593 def get_indices(self):
594 return np.unique(self.triples.ravel())
596 def todict(self):
597 return {'name': 'FixLinearTriatomic',
598 'kwargs': {'triples': self.triples.tolist()}}
600 def index_shuffle(self, atoms, ind):
601 """Shuffle the indices of the three atoms in this constraint"""
602 map = np.zeros(len(atoms), int)
603 map[ind] = 1
604 n = map.sum()
605 map[:] = -1
606 map[ind] = range(n)
607 triples = map[self.triples]
608 self.triples = triples[(triples != -1).all(1)]
609 if len(self.triples) == 0:
610 raise IndexError('Constraint not part of slice')
613class FixedMode(FixConstraint):
614 """Constrain atoms to move along directions orthogonal to
615 a given mode only. Initialize with a mode, such as one produced by
616 ase.vibrations.Vibrations.get_mode()."""
618 def __init__(self, mode):
619 mode = np.asarray(mode)
620 self.mode = (mode / np.sqrt((mode**2).sum())).reshape(-1)
622 def get_removed_dof(self, atoms):
623 return len(atoms)
625 def adjust_positions(self, atoms, newpositions):
626 newpositions = newpositions.ravel()
627 oldpositions = atoms.positions.ravel()
628 step = newpositions - oldpositions
629 newpositions -= self.mode * np.dot(step, self.mode)
631 def adjust_forces(self, atoms, forces):
632 forces = forces.ravel()
633 forces -= self.mode * np.dot(forces, self.mode)
635 def index_shuffle(self, atoms, ind):
636 eps = 1e-12
637 mode = self.mode.reshape(-1, 3)
638 excluded = np.ones(len(mode), dtype=bool)
639 excluded[ind] = False
640 if (abs(mode[excluded]) > eps).any():
641 raise IndexError('All nonzero parts of mode not in slice')
642 self.mode = mode[ind].ravel()
644 def get_indices(self):
645 # This function will never properly work because it works on all
646 # atoms and it has no idea how to tell how many atoms it is
647 # attached to. If it is being used, surely the user knows
648 # everything is being constrained.
649 return []
651 def todict(self):
652 return {'name': 'FixedMode',
653 'kwargs': {'mode': self.mode.tolist()}}
655 def __repr__(self):
656 return f'FixedMode({self.mode.tolist()})'
659def _normalize(direction):
660 if np.shape(direction) != (3,):
661 raise ValueError("len(direction) is {len(direction)}. Has to be 3")
663 direction = np.asarray(direction) / np.linalg.norm(direction)
664 return direction
667class FixedPlane(IndexedConstraint):
668 """
669 Constraint object for fixing chosen atoms to only move in a plane.
671 The plane is defined by its normal vector *direction*
672 """
674 def __init__(self, indices, direction):
675 """Constrain chosen atoms.
677 Parameters
678 ----------
679 indices : int or list of int
680 Index or indices for atoms that should be constrained
681 direction : list of 3 int
682 Direction of the normal vector
684 Examples
685 --------
686 Fix all Copper atoms to only move in the yz-plane:
688 >>> from ase.build import bulk
689 >>> from ase.constraints import FixedPlane
691 >>> atoms = bulk('Cu', 'fcc', a=3.6)
692 >>> c = FixedPlane(
693 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
694 ... direction=[1, 0, 0],
695 ... )
696 >>> atoms.set_constraint(c)
698 or constrain a single atom with the index 0 to move in the xy-plane:
700 >>> c = FixedPlane(indices=0, direction=[0, 0, 1])
701 >>> atoms.set_constraint(c)
702 """
703 super().__init__(indices=indices)
704 self.dir = _normalize(direction)
706 def adjust_positions(self, atoms, newpositions):
707 step = newpositions[self.index] - atoms.positions[self.index]
708 newpositions[self.index] -= _projection(step, self.dir)
710 def adjust_forces(self, atoms, forces):
711 forces[self.index] -= _projection(forces[self.index], self.dir)
713 def get_removed_dof(self, atoms):
714 return len(self.index)
716 def todict(self):
717 return {
718 'name': 'FixedPlane',
719 'kwargs': {'indices': self.index.tolist(),
720 'direction': self.dir.tolist()}
721 }
723 def __repr__(self):
724 return f'FixedPlane(indices={self.index}, {self.dir.tolist()})'
727def _projection(vectors, direction):
728 dotprods = vectors @ direction
729 projection = direction[None, :] * dotprods[:, None]
730 return projection
733class FixedLine(IndexedConstraint):
734 """
735 Constrain an atom index or a list of atom indices to move on a line only.
737 The line is defined by its vector *direction*
738 """
740 def __init__(self, indices, direction):
741 """Constrain chosen atoms.
743 Parameters
744 ----------
745 indices : int or list of int
746 Index or indices for atoms that should be constrained
747 direction : list of 3 int
748 Direction of the vector defining the line
750 Examples
751 --------
752 Fix all Copper atoms to only move in the x-direction:
754 >>> from ase.constraints import FixedLine
755 >>> c = FixedLine(
756 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
757 ... direction=[1, 0, 0],
758 ... )
759 >>> atoms.set_constraint(c)
761 or constrain a single atom with the index 0 to move in the z-direction:
763 >>> c = FixedLine(indices=0, direction=[0, 0, 1])
764 >>> atoms.set_constraint(c)
765 """
766 super().__init__(indices)
767 self.dir = _normalize(direction)
769 def adjust_positions(self, atoms, newpositions):
770 step = newpositions[self.index] - atoms.positions[self.index]
771 projection = _projection(step, self.dir)
772 newpositions[self.index] = atoms.positions[self.index] + projection
774 def adjust_forces(self, atoms, forces):
775 forces[self.index] = _projection(forces[self.index], self.dir)
777 def get_removed_dof(self, atoms):
778 return 2 * len(self.index)
780 def __repr__(self):
781 return f'FixedLine(indices={self.index}, {self.dir.tolist()})'
783 def todict(self):
784 return {
785 'name': 'FixedLine',
786 'kwargs': {'indices': self.index.tolist(),
787 'direction': self.dir.tolist()}
788 }
791class FixCartesian(IndexedConstraint):
792 'Fix an atom index *a* in the directions of the cartesian coordinates.'
794 def __init__(self, a, mask=(1, 1, 1)):
795 super().__init__(indices=a)
796 self.mask = ~np.asarray(mask, bool)
798 def get_removed_dof(self, atoms):
799 return (3 - self.mask.sum()) * len(self.index)
801 def adjust_positions(self, atoms, new):
802 step = new[self.index] - atoms.positions[self.index]
803 step *= self.mask[None, :]
804 new[self.index] = atoms.positions[self.index] + step
806 def adjust_forces(self, atoms, forces):
807 forces[self.index] *= self.mask[None, :]
809 def __repr__(self):
810 return 'FixCartesian(indices={}, mask={})'.format(
811 self.index.tolist(), list(~self.mask))
813 def todict(self):
814 return {'name': 'FixCartesian',
815 'kwargs': {'a': self.index.tolist(),
816 'mask': (~self.mask).tolist()}}
819class FixScaled(IndexedConstraint):
820 'Fix an atom index *a* in the directions of the unit vectors.'
822 def __init__(self, a, mask=(1, 1, 1), cell=None):
823 # XXX The unused cell keyword is there for compatibility
824 # with old trajectory files.
825 super().__init__(a)
826 self.mask = np.array(mask, bool)
828 def get_removed_dof(self, atoms):
829 return self.mask.sum() * len(self.index)
831 def adjust_positions(self, atoms, new):
832 cell = atoms.cell
833 scaled_old = cell.scaled_positions(atoms.positions[self.index])
834 scaled_new = cell.scaled_positions(new[self.index])
835 scaled_new[:, self.mask] = scaled_old[:, self.mask]
836 new[self.index] = cell.cartesian_positions(scaled_new)
838 def adjust_forces(self, atoms, forces):
839 # Forces are covariant to the coordinate transformation,
840 # use the inverse transformations
841 cell = atoms.cell
842 scaled_forces = cell.cartesian_positions(forces[self.index])
843 scaled_forces *= -(self.mask - 1)
844 forces[self.index] = cell.scaled_positions(scaled_forces)
846 def todict(self):
847 return {'name': 'FixScaled',
848 'kwargs': {'a': self.index.tolist(),
849 'mask': self.mask.tolist()}}
851 def __repr__(self):
852 return f'FixScaled({self.index.tolist()}, {self.mask})'
855# TODO: Better interface might be to use dictionaries in place of very
856# nested lists/tuples
857class FixInternals(FixConstraint):
858 """Constraint object for fixing multiple internal coordinates.
860 Allows fixing bonds, angles, dihedrals as well as linear combinations
861 of bonds (bondcombos).
863 Please provide angular units in degrees using `angles_deg` and
864 `dihedrals_deg`.
865 Fixing planar angles is not supported at the moment.
866 """
868 def __init__(self, bonds=None, angles=None, dihedrals=None,
869 angles_deg=None, dihedrals_deg=None,
870 bondcombos=None,
871 mic=False, epsilon=1.e-7):
872 """
873 A constrained internal coordinate is defined as a nested list:
874 '[value, [atom indices]]'. The constraint is initialized with a list of
875 constrained internal coordinates, i.e. '[[value, [atom indices]], ...]'.
876 If 'value' is None, the current value of the coordinate is constrained.
878 Parameters
879 ----------
880 bonds: nested python list, optional
881 List with targetvalue and atom indices defining the fixed bonds,
882 i.e. [[targetvalue, [index0, index1]], ...]
884 angles_deg: nested python list, optional
885 List with targetvalue and atom indices defining the fixedangles,
886 i.e. [[targetvalue, [index0, index1, index3]], ...]
888 dihedrals_deg: nested python list, optional
889 List with targetvalue and atom indices defining the fixed dihedrals,
890 i.e. [[targetvalue, [index0, index1, index3]], ...]
892 bondcombos: nested python list, optional
893 List with targetvalue, atom indices and linear coefficient defining
894 the fixed linear combination of bonds,
895 i.e. [[targetvalue, [[index0, index1, coefficient_for_bond],
896 [index1, index2, coefficient_for_bond]]], ...]
898 mic: bool, optional, default: False
899 Minimum image convention.
901 epsilon: float, optional, default: 1e-7
902 Convergence criterion.
903 """
904 warn_msg = 'Please specify {} in degrees using the {} argument.'
905 if angles:
906 warn(warn_msg.format('angles', 'angle_deg'), FutureWarning)
907 angles = np.asarray(angles)
908 angles[:, 0] = angles[:, 0] / np.pi * 180
909 angles = angles.tolist()
910 else:
911 angles = angles_deg
912 if dihedrals:
913 warn(warn_msg.format('dihedrals', 'dihedrals_deg'), FutureWarning)
914 dihedrals = np.asarray(dihedrals)
915 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180
916 dihedrals = dihedrals.tolist()
917 else:
918 dihedrals = dihedrals_deg
920 self.bonds = bonds or []
921 self.angles = angles or []
922 self.dihedrals = dihedrals or []
923 self.bondcombos = bondcombos or []
924 self.mic = mic
925 self.epsilon = epsilon
927 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)
928 + len(self.bondcombos))
930 # Initialize these at run-time:
931 self.constraints = []
932 self.initialized = False
934 def get_removed_dof(self, atoms):
935 return self.n
937 def initialize(self, atoms):
938 if self.initialized:
939 return
940 masses = np.repeat(atoms.get_masses(), 3)
941 cell = None
942 pbc = None
943 if self.mic:
944 cell = atoms.cell
945 pbc = atoms.pbc
946 self.constraints = []
947 for data, ConstrClass in [(self.bonds, self.FixBondLengthAlt),
948 (self.angles, self.FixAngle),
949 (self.dihedrals, self.FixDihedral),
950 (self.bondcombos, self.FixBondCombo)]:
951 for datum in data:
952 targetvalue = datum[0]
953 if targetvalue is None: # set to current value
954 targetvalue = ConstrClass.get_value(atoms, datum[1],
955 self.mic)
956 constr = ConstrClass(targetvalue, datum[1], masses, cell, pbc)
957 self.constraints.append(constr)
958 self.initialized = True
960 @staticmethod
961 def get_bondcombo(atoms, indices, mic=False):
962 """Convenience function to return the value of the bondcombo coordinate
963 (linear combination of bond lengths) for the given Atoms object 'atoms'.
964 Example: Get the current value of the linear combination of two bond
965 lengths defined as `bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]`."""
966 c = sum(df[2] * atoms.get_distance(*df[:2], mic=mic) for df in indices)
967 return c
969 def get_subconstraint(self, atoms, definition):
970 """Get pointer to a specific subconstraint.
971 Identification by its definition via indices (and coefficients)."""
972 self.initialize(atoms)
973 for subconstr in self.constraints:
974 if isinstance(definition[0], Sequence): # Combo constraint
975 defin = [d + [c] for d, c in zip(subconstr.indices,
976 subconstr.coefs)]
977 if defin == definition:
978 return subconstr
979 else: # identify primitive constraints by their indices
980 if subconstr.indices == [definition]:
981 return subconstr
982 raise ValueError('Given `definition` not found on Atoms object.')
984 def shuffle_definitions(self, shuffle_dic, internal_type):
985 dfns = [] # definitions
986 for dfn in internal_type: # e.g. for bond in self.bonds
987 append = True
988 new_dfn = [dfn[0], list(dfn[1])]
989 for old in dfn[1]:
990 if old in shuffle_dic:
991 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]
992 else:
993 append = False
994 break
995 if append:
996 dfns.append(new_dfn)
997 return dfns
999 def shuffle_combos(self, shuffle_dic, internal_type):
1000 dfns = [] # definitions
1001 for dfn in internal_type: # i.e. for bondcombo in self.bondcombos
1002 append = True
1003 all_indices = [idx[0:-1] for idx in dfn[1]]
1004 new_dfn = [dfn[0], list(dfn[1])]
1005 for i, indices in enumerate(all_indices):
1006 for old in indices:
1007 if old in shuffle_dic:
1008 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]
1009 else:
1010 append = False
1011 break
1012 if not append:
1013 break
1014 if append:
1015 dfns.append(new_dfn)
1016 return dfns
1018 def index_shuffle(self, atoms, ind):
1019 # See docstring of superclass
1020 self.initialize(atoms)
1021 shuffle_dic = dict(slice2enlist(ind, len(atoms)))
1022 shuffle_dic = {old: new for new, old in shuffle_dic.items()}
1023 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)
1024 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)
1025 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)
1026 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)
1027 self.initialized = False
1028 self.initialize(atoms)
1029 if len(self.constraints) == 0:
1030 raise IndexError('Constraint not part of slice')
1032 def get_indices(self):
1033 cons = []
1034 for dfn in self.bonds + self.dihedrals + self.angles:
1035 cons.extend(dfn[1])
1036 for dfn in self.bondcombos:
1037 for partial_dfn in dfn[1]:
1038 cons.extend(partial_dfn[0:-1]) # last index is the coefficient
1039 return list(set(cons))
1041 def todict(self):
1042 return {'name': 'FixInternals',
1043 'kwargs': {'bonds': self.bonds,
1044 'angles_deg': self.angles,
1045 'dihedrals_deg': self.dihedrals,
1046 'bondcombos': self.bondcombos,
1047 'mic': self.mic,
1048 'epsilon': self.epsilon}}
1050 def adjust_positions(self, atoms, newpos):
1051 self.initialize(atoms)
1052 for constraint in self.constraints:
1053 constraint.setup_jacobian(atoms.positions)
1054 for j in range(50):
1055 maxerr = 0.0
1056 for constraint in self.constraints:
1057 constraint.adjust_positions(atoms.positions, newpos)
1058 maxerr = max(abs(constraint.sigma), maxerr)
1059 if maxerr < self.epsilon:
1060 return
1061 msg = 'FixInternals.adjust_positions did not converge.'
1062 if any(constr.targetvalue > 175. or constr.targetvalue < 5. for constr
1063 in self.constraints if isinstance(constr, self.FixAngle)):
1064 msg += (' This may be caused by an almost planar angle.'
1065 ' Support for planar angles would require the'
1066 ' implementation of ghost, i.e. dummy, atoms.'
1067 ' See issue #868.')
1068 raise ValueError(msg)
1070 def adjust_forces(self, atoms, forces):
1071 """Project out translations and rotations and all other constraints"""
1072 self.initialize(atoms)
1073 positions = atoms.positions
1074 N = len(forces)
1075 list2_constraints = list(np.zeros((6, N, 3)))
1076 tx, ty, tz, rx, ry, rz = list2_constraints
1078 list_constraints = [r.ravel() for r in list2_constraints]
1080 tx[:, 0] = 1.0
1081 ty[:, 1] = 1.0
1082 tz[:, 2] = 1.0
1083 ff = forces.ravel()
1085 # Calculate the center of mass
1086 center = positions.sum(axis=0) / N
1088 rx[:, 1] = -(positions[:, 2] - center[2])
1089 rx[:, 2] = positions[:, 1] - center[1]
1090 ry[:, 0] = positions[:, 2] - center[2]
1091 ry[:, 2] = -(positions[:, 0] - center[0])
1092 rz[:, 0] = -(positions[:, 1] - center[1])
1093 rz[:, 1] = positions[:, 0] - center[0]
1095 # Normalizing transl., rotat. constraints
1096 for r in list2_constraints:
1097 r /= np.linalg.norm(r.ravel())
1099 # Add all angle, etc. constraint vectors
1100 for constraint in self.constraints:
1101 constraint.setup_jacobian(positions)
1102 constraint.adjust_forces(positions, forces)
1103 list_constraints.insert(0, constraint.jacobian)
1104 # QR DECOMPOSITION - GRAM SCHMIDT
1106 list_constraints = [r.ravel() for r in list_constraints]
1107 aa = np.column_stack(list_constraints)
1108 (aa, bb) = np.linalg.qr(aa)
1109 # Projection
1110 hh = []
1111 for i, constraint in enumerate(self.constraints):
1112 hh.append(aa[:, i] * np.row_stack(aa[:, i]))
1114 txx = aa[:, self.n] * np.row_stack(aa[:, self.n])
1115 tyy = aa[:, self.n + 1] * np.row_stack(aa[:, self.n + 1])
1116 tzz = aa[:, self.n + 2] * np.row_stack(aa[:, self.n + 2])
1117 rxx = aa[:, self.n + 3] * np.row_stack(aa[:, self.n + 3])
1118 ryy = aa[:, self.n + 4] * np.row_stack(aa[:, self.n + 4])
1119 rzz = aa[:, self.n + 5] * np.row_stack(aa[:, self.n + 5])
1120 T = txx + tyy + tzz + rxx + ryy + rzz
1121 for vec in hh:
1122 T += vec
1123 ff = np.dot(T, np.row_stack(ff))
1124 forces[:, :] -= np.dot(T, np.row_stack(ff)).reshape(-1, 3)
1126 def __repr__(self):
1127 constraints = [repr(constr) for constr in self.constraints]
1128 return f'FixInternals(_copy_init={constraints}, epsilon={self.epsilon})'
1130 # Classes for internal use in FixInternals
1131 class FixInternalsBase:
1132 """Base class for subclasses of FixInternals."""
1134 def __init__(self, targetvalue, indices, masses, cell, pbc):
1135 self.targetvalue = targetvalue # constant target value
1136 self.indices = [defin[0:-1] for defin in indices] # indices, defs
1137 self.coefs = np.asarray([defin[-1] for defin in indices])
1138 self.masses = masses
1139 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix
1140 self.sigma = 1. # difference between current and target value
1141 self.projected_force = None # helps optimizers scan along constr.
1142 self.cell = cell
1143 self.pbc = pbc
1145 def finalize_jacobian(self, pos, n_internals, n, derivs):
1146 """Populate jacobian with derivatives for `n_internals` defined
1147 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""
1148 jacobian = np.zeros((n_internals, *pos.shape))
1149 for i, idx in enumerate(self.indices):
1150 for j in range(n):
1151 jacobian[i, idx[j]] = derivs[i, j]
1152 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))
1153 return self.coefs @ jacobian
1155 def finalize_positions(self, newpos):
1156 jacobian = self.jacobian / self.masses
1157 lamda = -self.sigma / (jacobian @ self.get_jacobian(newpos))
1158 dnewpos = lamda * jacobian
1159 newpos += dnewpos.reshape(newpos.shape)
1161 def adjust_forces(self, positions, forces):
1162 self.projected_forces = ((self.jacobian @ forces.ravel())
1163 * self.jacobian)
1164 self.jacobian /= np.linalg.norm(self.jacobian)
1166 class FixBondCombo(FixInternalsBase):
1167 """Constraint subobject for fixing linear combination of bond lengths
1168 within FixInternals.
1170 sum_i( coef_i * bond_length_i ) = constant
1171 """
1173 def get_jacobian(self, pos):
1174 bondvectors = [pos[k] - pos[h] for h, k in self.indices]
1175 derivs = get_distances_derivatives(bondvectors, cell=self.cell,
1176 pbc=self.pbc)
1177 return self.finalize_jacobian(pos, len(bondvectors), 2, derivs)
1179 def setup_jacobian(self, pos):
1180 self.jacobian = self.get_jacobian(pos)
1182 def adjust_positions(self, oldpos, newpos):
1183 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]
1184 (_, ), (dists, ) = conditional_find_mic([bondvectors],
1185 cell=self.cell,
1186 pbc=self.pbc)
1187 value = self.coefs @ dists
1188 self.sigma = value - self.targetvalue
1189 self.finalize_positions(newpos)
1191 @staticmethod
1192 def get_value(atoms, indices, mic):
1193 return FixInternals.get_bondcombo(atoms, indices, mic)
1195 def __repr__(self):
1196 return (f'FixBondCombo({self.targetvalue}, {self.indices}, '
1197 '{self.coefs})')
1199 class FixBondLengthAlt(FixBondCombo):
1200 """Constraint subobject for fixing bond length within FixInternals.
1201 Fix distance between atoms with indices a1, a2."""
1203 def __init__(self, targetvalue, indices, masses, cell, pbc):
1204 if targetvalue <= 0.:
1205 raise ZeroDivisionError('Invalid targetvalue for fixed bond')
1206 indices = [list(indices) + [1.]] # bond definition with coef 1.
1207 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1209 @staticmethod
1210 def get_value(atoms, indices, mic):
1211 return atoms.get_distance(*indices, mic=mic)
1213 def __repr__(self):
1214 return f'FixBondLengthAlt({self.targetvalue}, {self.indices})'
1216 class FixAngle(FixInternalsBase):
1217 """Constraint subobject for fixing an angle within FixInternals.
1219 Convergence is potentially problematic for angles very close to
1220 0 or 180 degrees as there is a singularity in the Cartesian derivative.
1221 Fixing planar angles is therefore not supported at the moment.
1222 """
1224 def __init__(self, targetvalue, indices, masses, cell, pbc):
1225 """Fix atom movement to construct a constant angle."""
1226 if targetvalue <= 0. or targetvalue >= 180.:
1227 raise ZeroDivisionError('Invalid targetvalue for fixed angle')
1228 indices = [list(indices) + [1.]] # angle definition with coef 1.
1229 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1231 def gather_vectors(self, pos):
1232 v0 = [pos[h] - pos[k] for h, k, l in self.indices]
1233 v1 = [pos[l] - pos[k] for h, k, l in self.indices]
1234 return v0, v1
1236 def get_jacobian(self, pos):
1237 v0, v1 = self.gather_vectors(pos)
1238 derivs = get_angles_derivatives(v0, v1, cell=self.cell,
1239 pbc=self.pbc)
1240 return self.finalize_jacobian(pos, len(v0), 3, derivs)
1242 def setup_jacobian(self, pos):
1243 self.jacobian = self.get_jacobian(pos)
1245 def adjust_positions(self, oldpos, newpos):
1246 v0, v1 = self.gather_vectors(newpos)
1247 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)
1248 self.sigma = value - self.targetvalue
1249 self.finalize_positions(newpos)
1251 @staticmethod
1252 def get_value(atoms, indices, mic):
1253 return atoms.get_angle(*indices, mic=mic)
1255 def __repr__(self):
1256 return f'FixAngle({self.targetvalue}, {self.indices})'
1258 class FixDihedral(FixInternalsBase):
1259 """Constraint subobject for fixing a dihedral angle within FixInternals.
1261 A dihedral becomes undefined when at least one of the inner two angles
1262 becomes planar. Make sure to avoid this situation.
1263 """
1265 def __init__(self, targetvalue, indices, masses, cell, pbc):
1266 indices = [list(indices) + [1.]] # dihedral def. with coef 1.
1267 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1269 def gather_vectors(self, pos):
1270 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]
1271 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]
1272 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]
1273 return v0, v1, v2
1275 def get_jacobian(self, pos):
1276 v0, v1, v2 = self.gather_vectors(pos)
1277 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,
1278 pbc=self.pbc)
1279 return self.finalize_jacobian(pos, len(v0), 4, derivs)
1281 def setup_jacobian(self, pos):
1282 self.jacobian = self.get_jacobian(pos)
1284 def adjust_positions(self, oldpos, newpos):
1285 v0, v1, v2 = self.gather_vectors(newpos)
1286 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1287 # apply minimum dihedral difference 'convention': (diff <= 180)
1288 self.sigma = (value - self.targetvalue + 180) % 360 - 180
1289 self.finalize_positions(newpos)
1291 @staticmethod
1292 def get_value(atoms, indices, mic):
1293 return atoms.get_dihedral(*indices, mic=mic)
1295 def __repr__(self):
1296 return f'FixDihedral({self.targetvalue}, {self.indices})'
1299class FixParametricRelations(FixConstraint):
1301 def __init__(
1302 self,
1303 indices,
1304 Jacobian,
1305 const_shift,
1306 params=None,
1307 eps=1e-12,
1308 use_cell=False,
1309 ):
1310 """Constrains the degrees of freedom to act in a reduced parameter
1311 space defined by the Jacobian
1313 These constraints are based off the work in:
1314 https://arxiv.org/abs/1908.01610
1316 The constraints linearly maps the full 3N degrees of freedom,
1317 where N is number of active lattice vectors/atoms onto a
1318 reduced subset of M free parameters, where M <= 3*N. The
1319 Jacobian matrix and constant shift vector map the full set of
1320 degrees of freedom onto the reduced parameter space.
1322 Currently the constraint is set up to handle either atomic
1323 positions or lattice vectors at one time, but not both. To do
1324 both simply add a two constraints for each set. This is done
1325 to keep the mathematics behind the operations separate.
1327 It would be possible to extend these constraints to allow
1328 non-linear transformations if functionality to update the
1329 Jacobian at each position update was included. This would
1330 require passing an update function evaluate it every time
1331 adjust_positions is callled. This is currently NOT supported,
1332 and there are no plans to implement it in the future.
1334 Args:
1335 indices (list of int): indices of the constrained atoms
1336 (if not None or empty then cell_indices must be None or Empty)
1337 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))):
1338 The Jacobian describing
1339 the parameter space transformation
1340 const_shift (np.ndarray(shape=(3*len(indices)))):
1341 A vector describing the constant term
1342 in the transformation not accounted for in the Jacobian
1343 params (list of str):
1344 parameters used in the parametric representation
1345 if None a list is generated based on the shape of the Jacobian
1346 eps (float): a small number to compare the similarity of
1347 numbers and set the precision used
1348 to generate the constraint expressions
1349 use_cell (bool): if True then act on the cell object
1351 """
1352 self.indices = np.array(indices)
1353 self.Jacobian = np.array(Jacobian)
1354 self.const_shift = np.array(const_shift)
1356 assert self.const_shift.shape[0] == 3 * len(self.indices)
1357 assert self.Jacobian.shape[0] == 3 * len(self.indices)
1359 self.eps = eps
1360 self.use_cell = use_cell
1362 if params is None:
1363 params = []
1364 if self.Jacobian.shape[1] > 0:
1365 int_fmt_str = "{:0" + \
1366 str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"
1367 for param_ind in range(self.Jacobian.shape[1]):
1368 params.append("param_" + int_fmt_str.format(param_ind))
1369 else:
1370 assert len(params) == self.Jacobian.shape[-1]
1372 self.params = params
1374 self.Jacobian_inv = np.linalg.inv(
1375 self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T
1377 @classmethod
1378 def from_expressions(cls, indices, params, expressions,
1379 eps=1e-12, use_cell=False):
1380 """Converts the expressions into a Jacobian Matrix/const_shift
1381 vector and constructs a FixParametricRelations constraint
1383 The expressions must be a list like object of size 3*N and
1384 elements must be ordered as:
1385 [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],
1386 where i, j, and k are the first, second and third
1387 component of the atomic position/lattice
1388 vector. Currently only linear operations are allowed to be
1389 included in the expressions so
1390 only terms like:
1391 - const * param_0
1392 - sqrt[const] * param_1
1393 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M
1394 where const is any real number and param_0, param_1, ..., param_M are
1395 the parameters passed in
1396 params, are allowed.
1398 For example, fractional atomic position constraints for wurtzite are:
1399 params = ["z1", "z2"]
1400 expressions = [
1401 "1.0/3.0", "2.0/3.0", "z1",
1402 "2.0/3.0", "1.0/3.0", "0.5 + z1",
1403 "1.0/3.0", "2.0/3.0", "z2",
1404 "2.0/3.0", "1.0/3.0", "0.5 + z2",
1405 ]
1407 For diamond are:
1408 params = []
1409 expressions = [
1410 "0.0", "0.0", "0.0",
1411 "0.25", "0.25", "0.25",
1412 ],
1414 and for stannite are
1415 params=["x4", "z4"]
1416 expressions = [
1417 "0.0", "0.0", "0.0",
1418 "0.0", "0.5", "0.5",
1419 "0.75", "0.25", "0.5",
1420 "0.25", "0.75", "0.5",
1421 "x4 + z4", "x4 + z4", "2*x4",
1422 "x4 - z4", "x4 - z4", "-2*x4",
1423 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",
1424 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",
1425 ]
1427 Args:
1428 indices (list of int): indices of the constrained atoms
1429 (if not None or empty then cell_indices must be None or Empty)
1430 params (list of str): parameters used in the
1431 parametric representation
1432 expressions (list of str): expressions used to convert from the
1433 parametric to the real space representation
1434 eps (float): a small number to compare the similarity of
1435 numbers and set the precision used
1436 to generate the constraint expressions
1437 use_cell (bool): if True then act on the cell object
1439 Returns:
1440 cls(
1441 indices,
1442 Jacobian generated from expressions,
1443 const_shift generated from expressions,
1444 params,
1445 eps-12,
1446 use_cell,
1447 )
1448 """
1449 Jacobian = np.zeros((3 * len(indices), len(params)))
1450 const_shift = np.zeros(3 * len(indices))
1452 for expr_ind, expression in enumerate(expressions):
1453 expression = expression.strip()
1455 # Convert subtraction to addition
1456 expression = expression.replace("-", "+(-1.0)*")
1457 if "+" == expression[0]:
1458 expression = expression[1:]
1459 elif "(+" == expression[:2]:
1460 expression = "(" + expression[2:]
1462 # Explicitly add leading zeros so when replacing param_1 with 0.0
1463 # param_11 does not become 0.01
1464 int_fmt_str = "{:0" + \
1465 str(int(np.ceil(np.log10(len(params) + 1)))) + "d}"
1467 param_dct = {}
1468 param_map = {}
1470 # Construct a standardized param template for A/B filling
1471 for param_ind, param in enumerate(params):
1472 param_str = "param_" + int_fmt_str.format(param_ind)
1473 param_map[param] = param_str
1474 param_dct[param_str] = 0.0
1476 # Replace the parameters according to the map
1477 # Sort by string length (long to short) to prevent cases like x11
1478 # becoming f"{param_map["x1"]}1"
1479 for param in sorted(params, key=lambda s: -1.0 * len(s)):
1480 expression = expression.replace(param, param_map[param])
1482 # Partial linearity check
1483 for express_sec in expression.split("+"):
1484 in_sec = [param in express_sec for param in param_dct]
1485 n_params_in_sec = len(np.where(np.array(in_sec))[0])
1486 if n_params_in_sec > 1:
1487 raise ValueError(
1488 "FixParametricRelations expressions must be linear.")
1490 const_shift[expr_ind] = float(
1491 eval_expression(expression, param_dct))
1493 for param_ind in range(len(params)):
1494 param_str = "param_" + int_fmt_str.format(param_ind)
1495 if param_str not in expression:
1496 Jacobian[expr_ind, param_ind] = 0.0
1497 continue
1498 param_dct[param_str] = 1.0
1499 test_1 = float(eval_expression(expression, param_dct))
1500 test_1 -= const_shift[expr_ind]
1501 Jacobian[expr_ind, param_ind] = test_1
1503 param_dct[param_str] = 2.0
1504 test_2 = float(eval_expression(expression, param_dct))
1505 test_2 -= const_shift[expr_ind]
1506 if abs(test_2 / test_1 - 2.0) > eps:
1507 raise ValueError(
1508 "FixParametricRelations expressions must be linear.")
1509 param_dct[param_str] = 0.0
1511 args = [
1512 indices,
1513 Jacobian,
1514 const_shift,
1515 params,
1516 eps,
1517 use_cell,
1518 ]
1519 if cls is FixScaledParametricRelations:
1520 args = args[:-1]
1521 return cls(*args)
1523 @property
1524 def expressions(self):
1525 """Generate the expressions represented by the current self.Jacobian
1526 and self.const_shift objects"""
1527 expressions = []
1528 per = int(round(-1 * np.log10(self.eps)))
1529 fmt_str = "{:." + str(per + 1) + "g}"
1530 for index, shift_val in enumerate(self.const_shift):
1531 exp = ""
1532 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(
1533 shift_val) > self.eps:
1534 exp += fmt_str.format(shift_val)
1536 param_exp = ""
1537 for param_index, jacob_val in enumerate(self.Jacobian[index]):
1538 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)
1539 if abs_jacob_val < self.eps:
1540 continue
1542 param = self.params[param_index]
1543 if param_exp or exp:
1544 if jacob_val > -1.0 * self.eps:
1545 param_exp += " + "
1546 else:
1547 param_exp += " - "
1548 elif (not exp) and (not param_exp) and (
1549 jacob_val < -1.0 * self.eps):
1550 param_exp += "-"
1552 if np.abs(abs_jacob_val - 1.0) <= self.eps:
1553 param_exp += f"{param:s}"
1554 else:
1555 param_exp += (fmt_str +
1556 "*{:s}").format(abs_jacob_val, param)
1558 exp += param_exp
1560 expressions.append(exp)
1561 return np.array(expressions).reshape((-1, 3))
1563 def todict(self):
1564 """Create a dictionary representation of the constraint"""
1565 return {
1566 "name": type(self).__name__,
1567 "kwargs": {
1568 "indices": self.indices,
1569 "params": self.params,
1570 "Jacobian": self.Jacobian,
1571 "const_shift": self.const_shift,
1572 "eps": self.eps,
1573 "use_cell": self.use_cell,
1574 }
1575 }
1577 def __repr__(self):
1578 """The str representation of the constraint"""
1579 if len(self.indices) > 1:
1580 indices_str = "[{:d}, ..., {:d}]".format(
1581 self.indices[0], self.indices[-1])
1582 else:
1583 indices_str = f"[{self.indices[0]:d}]"
1585 if len(self.params) > 1:
1586 params_str = "[{:s}, ..., {:s}]".format(
1587 self.params[0], self.params[-1])
1588 elif len(self.params) == 1:
1589 params_str = f"[{self.params[0]:s}]"
1590 else:
1591 params_str = "[]"
1593 return '{:s}({:s}, {:s}, ..., {:e})'.format(
1594 type(self).__name__,
1595 indices_str,
1596 params_str,
1597 self.eps
1598 )
1601class FixScaledParametricRelations(FixParametricRelations):
1603 def __init__(
1604 self,
1605 indices,
1606 Jacobian,
1607 const_shift,
1608 params=None,
1609 eps=1e-12,
1610 ):
1611 """The fractional coordinate version of FixParametricRelations
1613 All arguments are the same, but since this is for fractional
1614 coordinates use_cell is false"""
1615 super().__init__(
1616 indices,
1617 Jacobian,
1618 const_shift,
1619 params,
1620 eps,
1621 False,
1622 )
1624 def adjust_contravariant(self, cell, vecs, B):
1625 """Adjust the values of a set of vectors that are contravariant
1626 with the unit transformation"""
1627 scaled = cell.scaled_positions(vecs).flatten()
1628 scaled = self.Jacobian_inv @ (scaled - B)
1629 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))
1631 return cell.cartesian_positions(scaled)
1633 def adjust_positions(self, atoms, positions):
1634 """Adjust positions of the atoms to match the constraints"""
1635 positions[self.indices] = self.adjust_contravariant(
1636 atoms.cell,
1637 positions[self.indices],
1638 self.const_shift,
1639 )
1640 positions[self.indices] = self.adjust_B(
1641 atoms.cell, positions[self.indices])
1643 def adjust_B(self, cell, positions):
1644 """Wraps the positions back to the unit cell and adjust B to
1645 keep track of this change"""
1646 fractional = cell.scaled_positions(positions)
1647 wrapped_fractional = (fractional % 1.0) % 1.0
1648 self.const_shift += np.round(wrapped_fractional - fractional).flatten()
1649 return cell.cartesian_positions(wrapped_fractional)
1651 def adjust_momenta(self, atoms, momenta):
1652 """Adjust momenta of the atoms to match the constraints"""
1653 momenta[self.indices] = self.adjust_contravariant(
1654 atoms.cell,
1655 momenta[self.indices],
1656 np.zeros(self.const_shift.shape),
1657 )
1659 def adjust_forces(self, atoms, forces):
1660 """Adjust forces of the atoms to match the constraints"""
1661 # Forces are coavarient to the coordinate transformation, use the
1662 # inverse transformations
1663 cart2frac_jacob = np.zeros(2 * (3 * len(atoms),))
1664 for i_atom in range(len(atoms)):
1665 cart2frac_jacob[3 * i_atom:3 * (i_atom + 1),
1666 3 * i_atom:3 * (i_atom + 1)] = atoms.cell.T
1668 jacobian = cart2frac_jacob @ self.Jacobian
1669 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T
1671 reduced_forces = jacobian.T @ forces.flatten()
1672 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)
1674 def todict(self):
1675 """Create a dictionary representation of the constraint"""
1676 dct = super().todict()
1677 del dct["kwargs"]["use_cell"]
1678 return dct
1681class FixCartesianParametricRelations(FixParametricRelations):
1683 def __init__(
1684 self,
1685 indices,
1686 Jacobian,
1687 const_shift,
1688 params=None,
1689 eps=1e-12,
1690 use_cell=False,
1691 ):
1692 """The Cartesian coordinate version of FixParametricRelations"""
1693 super().__init__(
1694 indices,
1695 Jacobian,
1696 const_shift,
1697 params,
1698 eps,
1699 use_cell,
1700 )
1702 def adjust_contravariant(self, vecs, B):
1703 """Adjust the values of a set of vectors that are contravariant with
1704 the unit transformation"""
1705 vecs = self.Jacobian_inv @ (vecs.flatten() - B)
1706 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))
1707 return vecs
1709 def adjust_positions(self, atoms, positions):
1710 """Adjust positions of the atoms to match the constraints"""
1711 if self.use_cell:
1712 return
1713 positions[self.indices] = self.adjust_contravariant(
1714 positions[self.indices],
1715 self.const_shift,
1716 )
1718 def adjust_momenta(self, atoms, momenta):
1719 """Adjust momenta of the atoms to match the constraints"""
1720 if self.use_cell:
1721 return
1722 momenta[self.indices] = self.adjust_contravariant(
1723 momenta[self.indices],
1724 np.zeros(self.const_shift.shape),
1725 )
1727 def adjust_forces(self, atoms, forces):
1728 """Adjust forces of the atoms to match the constraints"""
1729 if self.use_cell:
1730 return
1732 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()
1733 forces[self.indices] = (self.Jacobian_inv.T @
1734 forces_reduced).reshape(-1, 3)
1736 def adjust_cell(self, atoms, cell):
1737 """Adjust the cell of the atoms to match the constraints"""
1738 if not self.use_cell:
1739 return
1740 cell[self.indices] = self.adjust_contravariant(
1741 cell[self.indices],
1742 np.zeros(self.const_shift.shape),
1743 )
1745 def adjust_stress(self, atoms, stress):
1746 """Adjust the stress of the atoms to match the constraints"""
1747 if not self.use_cell:
1748 return
1750 stress_3x3 = voigt_6_to_full_3x3_stress(stress)
1751 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()
1752 stress_3x3[self.indices] = (
1753 self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)
1755 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)
1758class Hookean(FixConstraint):
1759 """Applies a Hookean restorative force between a pair of atoms, an atom
1760 and a point, or an atom and a plane."""
1762 def __init__(self, a1, a2, k, rt=None):
1763 """Forces two atoms to stay close together by applying no force if
1764 they are below a threshold length, rt, and applying a Hookean
1765 restorative force when the distance between them exceeds rt. Can
1766 also be used to tether an atom to a fixed point in space or to a
1767 distance above a plane.
1769 a1 : int
1770 Index of atom 1
1771 a2 : one of three options
1772 1) index of atom 2
1773 2) a fixed point in cartesian space to which to tether a1
1774 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.
1775 k : float
1776 Hooke's law (spring) constant to apply when distance
1777 exceeds threshold_length. Units of eV A^-2.
1778 rt : float
1779 The threshold length below which there is no force. The
1780 length is 1) between two atoms, 2) between atom and point.
1781 This argument is not supplied in case 3. Units of A.
1783 If a plane is specified, the Hooke's law force is applied if the atom
1784 is on the normal side of the plane. For instance, the plane with
1785 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z
1786 intercept of +7 and a normal vector pointing in the +z direction.
1787 If the atom has z > 7, then a downward force would be applied of
1788 k * (atom.z - 7). The same plane with the normal vector pointing in
1789 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).
1791 References:
1793 Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014)
1794 https://link.springer.com/article/10.1007%2Fs11244-013-0161-8
1795 """
1797 if isinstance(a2, int):
1798 self._type = 'two atoms'
1799 self.indices = [a1, a2]
1800 elif len(a2) == 3:
1801 self._type = 'point'
1802 self.index = a1
1803 self.origin = np.array(a2)
1804 elif len(a2) == 4:
1805 self._type = 'plane'
1806 self.index = a1
1807 self.plane = a2
1808 else:
1809 raise RuntimeError('Unknown type for a2')
1810 self.threshold = rt
1811 self.spring = k
1813 def get_removed_dof(self, atoms):
1814 return 0
1816 def todict(self):
1817 dct = {'name': 'Hookean'}
1818 dct['kwargs'] = {'rt': self.threshold,
1819 'k': self.spring}
1820 if self._type == 'two atoms':
1821 dct['kwargs']['a1'] = self.indices[0]
1822 dct['kwargs']['a2'] = self.indices[1]
1823 elif self._type == 'point':
1824 dct['kwargs']['a1'] = self.index
1825 dct['kwargs']['a2'] = self.origin
1826 elif self._type == 'plane':
1827 dct['kwargs']['a1'] = self.index
1828 dct['kwargs']['a2'] = self.plane
1829 else:
1830 raise NotImplementedError(f'Bad type: {self._type}')
1831 return dct
1833 def adjust_positions(self, atoms, newpositions):
1834 pass
1836 def adjust_momenta(self, atoms, momenta):
1837 pass
1839 def adjust_forces(self, atoms, forces):
1840 positions = atoms.positions
1841 if self._type == 'plane':
1842 A, B, C, D = self.plane
1843 x, y, z = positions[self.index]
1844 d = ((A * x + B * y + C * z + D) /
1845 np.sqrt(A**2 + B**2 + C**2))
1846 if d < 0:
1847 return
1848 magnitude = self.spring * d
1849 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))
1850 forces[self.index] += direction * magnitude
1851 return
1852 if self._type == 'two atoms':
1853 p1, p2 = positions[self.indices]
1854 elif self._type == 'point':
1855 p1 = positions[self.index]
1856 p2 = self.origin
1857 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1858 bondlength = np.linalg.norm(displace)
1859 if bondlength > self.threshold:
1860 magnitude = self.spring * (bondlength - self.threshold)
1861 direction = displace / np.linalg.norm(displace)
1862 if self._type == 'two atoms':
1863 forces[self.indices[0]] += direction * magnitude
1864 forces[self.indices[1]] -= direction * magnitude
1865 else:
1866 forces[self.index] += direction * magnitude
1868 def adjust_potential_energy(self, atoms):
1869 """Returns the difference to the potential energy due to an active
1870 constraint. (That is, the quantity returned is to be added to the
1871 potential energy.)"""
1872 positions = atoms.positions
1873 if self._type == 'plane':
1874 A, B, C, D = self.plane
1875 x, y, z = positions[self.index]
1876 d = ((A * x + B * y + C * z + D) /
1877 np.sqrt(A**2 + B**2 + C**2))
1878 if d > 0:
1879 return 0.5 * self.spring * d**2
1880 else:
1881 return 0.
1882 if self._type == 'two atoms':
1883 p1, p2 = positions[self.indices]
1884 elif self._type == 'point':
1885 p1 = positions[self.index]
1886 p2 = self.origin
1887 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1888 bondlength = np.linalg.norm(displace)
1889 if bondlength > self.threshold:
1890 return 0.5 * self.spring * (bondlength - self.threshold)**2
1891 else:
1892 return 0.
1894 def get_indices(self):
1895 if self._type == 'two atoms':
1896 return self.indices
1897 elif self._type == 'point':
1898 return self.index
1899 elif self._type == 'plane':
1900 return self.index
1902 def index_shuffle(self, atoms, ind):
1903 # See docstring of superclass
1904 if self._type == 'two atoms':
1905 newa = [-1, -1] # Signal error
1906 for new, old in slice2enlist(ind, len(atoms)):
1907 for i, a in enumerate(self.indices):
1908 if old == a:
1909 newa[i] = new
1910 if newa[0] == -1 or newa[1] == -1:
1911 raise IndexError('Constraint not part of slice')
1912 self.indices = newa
1913 elif (self._type == 'point') or (self._type == 'plane'):
1914 newa = -1 # Signal error
1915 for new, old in slice2enlist(ind, len(atoms)):
1916 if old == self.index:
1917 newa = new
1918 break
1919 if newa == -1:
1920 raise IndexError('Constraint not part of slice')
1921 self.index = newa
1923 def __repr__(self):
1924 if self._type == 'two atoms':
1925 return 'Hookean(%d, %d)' % tuple(self.indices)
1926 elif self._type == 'point':
1927 return 'Hookean(%d) to cartesian' % self.index
1928 else:
1929 return 'Hookean(%d) to plane' % self.index
1932class ExternalForce(FixConstraint):
1933 """Constraint object for pulling two atoms apart by an external force.
1935 You can combine this constraint for example with FixBondLength but make
1936 sure that *ExternalForce* comes first in the list if there are overlaps
1937 between atom1-2 and atom3-4:
1939 >>> from ase.build import bulk
1941 >>> atoms = bulk('Cu', 'fcc', a=3.6)
1942 >>> atom1, atom2, atom3, atom4 = atoms[:4]
1943 >>> fext = 1.0
1944 >>> con1 = ExternalForce(atom1, atom2, f_ext)
1945 >>> con2 = FixBondLength(atom3, atom4)
1946 >>> atoms.set_constraint([con1, con2])
1948 see ase/test/external_force.py"""
1950 def __init__(self, a1, a2, f_ext):
1951 self.indices = [a1, a2]
1952 self.external_force = f_ext
1954 def get_removed_dof(self, atoms):
1955 return 0
1957 def adjust_positions(self, atoms, new):
1958 pass
1960 def adjust_forces(self, atoms, forces):
1961 dist = np.subtract.reduce(atoms.positions[self.indices])
1962 force = self.external_force * dist / np.linalg.norm(dist)
1963 forces[self.indices] += (force, -force)
1965 def adjust_potential_energy(self, atoms):
1966 dist = np.subtract.reduce(atoms.positions[self.indices])
1967 return -np.linalg.norm(dist) * self.external_force
1969 def index_shuffle(self, atoms, ind):
1970 """Shuffle the indices of the two atoms in this constraint"""
1971 newa = [-1, -1] # Signal error
1972 for new, old in slice2enlist(ind, len(atoms)):
1973 for i, a in enumerate(self.indices):
1974 if old == a:
1975 newa[i] = new
1976 if newa[0] == -1 or newa[1] == -1:
1977 raise IndexError('Constraint not part of slice')
1978 self.indices = newa
1980 def __repr__(self):
1981 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],
1982 self.indices[1],
1983 self.external_force)
1985 def todict(self):
1986 return {'name': 'ExternalForce',
1987 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
1988 'f_ext': self.external_force}}
1991class MirrorForce(FixConstraint):
1992 """Constraint object for mirroring the force between two atoms.
1994 This class is designed to find a transition state with the help of a
1995 single optimization. It can be used if the transition state belongs to a
1996 bond breaking reaction. First the given bond length will be fixed until
1997 all other degrees of freedom are optimized, then the forces of the two
1998 atoms will be mirrored to find the transition state. The mirror plane is
1999 perpendicular to the connecting line of the atoms. Transition states in
2000 dependence of the force can be obtained by stretching the molecule and
2001 fixing its total length with *FixBondLength* or by using *ExternalForce*
2002 during the optimization with *MirrorForce*.
2004 Parameters
2005 ----------
2006 a1: int
2007 First atom index.
2008 a2: int
2009 Second atom index.
2010 max_dist: float
2011 Upper limit of the bond length interval where the transition state
2012 can be found.
2013 min_dist: float
2014 Lower limit of the bond length interval where the transition state
2015 can be found.
2016 fmax: float
2017 Maximum force used for the optimization.
2019 Notes
2020 -----
2021 You can combine this constraint for example with FixBondLength but make
2022 sure that *MirrorForce* comes first in the list if there are overlaps
2023 between atom1-2 and atom3-4:
2025 >>> from ase.build import bulk
2027 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2028 >>> atom1, atom2, atom3, atom4 = atoms[:4]
2029 >>> con1 = MirrorForce(atom1, atom2)
2030 >>> con2 = FixBondLength(atom3, atom4)
2031 >>> atoms.set_constraint([con1, con2])
2033 """
2035 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):
2036 self.indices = [a1, a2]
2037 self.min_dist = min_dist
2038 self.max_dist = max_dist
2039 self.fmax = fmax
2041 def adjust_positions(self, atoms, new):
2042 pass
2044 def adjust_forces(self, atoms, forces):
2045 dist = np.subtract.reduce(atoms.positions[self.indices])
2046 d = np.linalg.norm(dist)
2047 if (d < self.min_dist) or (d > self.max_dist):
2048 # Stop structure optimization
2049 forces[:] *= 0
2050 return
2051 dist /= d
2052 df = np.subtract.reduce(forces[self.indices])
2053 f = df.dot(dist)
2054 con_saved = atoms.constraints
2055 try:
2056 con = [con for con in con_saved
2057 if not isinstance(con, MirrorForce)]
2058 atoms.set_constraint(con)
2059 forces_copy = atoms.get_forces()
2060 finally:
2061 atoms.set_constraint(con_saved)
2062 df1 = -1 / 2. * f * dist
2063 forces_copy[self.indices] += (df1, -df1)
2064 # Check if forces would be converged if the bond with mirrored forces
2065 # would also be fixed
2066 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2067 factor = 1.
2068 else:
2069 factor = 0.
2070 df1 = -(1 + factor) / 2. * f * dist
2071 forces[self.indices] += (df1, -df1)
2073 def index_shuffle(self, atoms, ind):
2074 """Shuffle the indices of the two atoms in this constraint
2076 """
2077 newa = [-1, -1] # Signal error
2078 for new, old in slice2enlist(ind, len(atoms)):
2079 for i, a in enumerate(self.indices):
2080 if old == a:
2081 newa[i] = new
2082 if newa[0] == -1 or newa[1] == -1:
2083 raise IndexError('Constraint not part of slice')
2084 self.indices = newa
2086 def __repr__(self):
2087 return 'MirrorForce(%d, %d, %f, %f, %f)' % (
2088 self.indices[0], self.indices[1], self.max_dist, self.min_dist,
2089 self.fmax)
2091 def todict(self):
2092 return {'name': 'MirrorForce',
2093 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2094 'max_dist': self.max_dist,
2095 'min_dist': self.min_dist, 'fmax': self.fmax}}
2098class MirrorTorque(FixConstraint):
2099 """Constraint object for mirroring the torque acting on a dihedral
2100 angle defined by four atoms.
2102 This class is designed to find a transition state with the help of a
2103 single optimization. It can be used if the transition state belongs to a
2104 cis-trans-isomerization with a change of dihedral angle. First the given
2105 dihedral angle will be fixed until all other degrees of freedom are
2106 optimized, then the torque acting on the dihedral angle will be mirrored
2107 to find the transition state. Transition states in
2108 dependence of the force can be obtained by stretching the molecule and
2109 fixing its total length with *FixBondLength* or by using *ExternalForce*
2110 during the optimization with *MirrorTorque*.
2112 This constraint can be used to find
2113 transition states of cis-trans-isomerization.
2115 a1 a4
2116 | |
2117 a2 __ a3
2119 Parameters
2120 ----------
2121 a1: int
2122 First atom index.
2123 a2: int
2124 Second atom index.
2125 a3: int
2126 Third atom index.
2127 a4: int
2128 Fourth atom index.
2129 max_angle: float
2130 Upper limit of the dihedral angle interval where the transition state
2131 can be found.
2132 min_angle: float
2133 Lower limit of the dihedral angle interval where the transition state
2134 can be found.
2135 fmax: float
2136 Maximum force used for the optimization.
2138 Notes
2139 -----
2140 You can combine this constraint for example with FixBondLength but make
2141 sure that *MirrorTorque* comes first in the list if there are overlaps
2142 between atom1-4 and atom5-6:
2144 >>> from ase.build import bulk
2146 >>> atoms = bulk('Cu', 'fcc', a=3.6)
2147 >>> atom1, atom2, atom3, atom4, atom5, atom6 = atoms[:6]
2148 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)
2149 >>> con2 = FixBondLength(atom5, atom6)
2150 >>> atoms.set_constraint([con1, con2])
2152 """
2154 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,
2155 fmax=0.1):
2156 self.indices = [a1, a2, a3, a4]
2157 self.min_angle = min_angle
2158 self.max_angle = max_angle
2159 self.fmax = fmax
2161 def adjust_positions(self, atoms, new):
2162 pass
2164 def adjust_forces(self, atoms, forces):
2165 angle = atoms.get_dihedral(self.indices[0], self.indices[1],
2166 self.indices[2], self.indices[3])
2167 angle *= np.pi / 180.
2168 if (angle < self.min_angle) or (angle > self.max_angle):
2169 # Stop structure optimization
2170 forces[:] *= 0
2171 return
2172 p = atoms.positions[self.indices]
2173 f = forces[self.indices]
2175 f0 = (f[1] + f[2]) / 2.
2176 ff = f - f0
2177 p0 = (p[2] + p[1]) / 2.
2178 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)
2179 fff = ff - np.cross(m0, p - p0)
2180 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \
2181 (p[1] - p0).dot(p[1] - p0)
2182 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \
2183 (p[2] - p0).dot(p[2] - p0)
2184 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \
2185 np.linalg.norm(p[1] - p0)
2186 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \
2187 np.linalg.norm(p[2] - p0)
2188 omegap = omegap1 + omegap2
2189 con_saved = atoms.constraints
2190 try:
2191 con = [con for con in con_saved
2192 if not isinstance(con, MirrorTorque)]
2193 atoms.set_constraint(con)
2194 forces_copy = atoms.get_forces()
2195 finally:
2196 atoms.set_constraint(con_saved)
2197 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \
2198 np.linalg.norm(p[1] - p0)
2199 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \
2200 np.linalg.norm(p[2] - p0)
2201 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2202 # Check if forces would be converged if the dihedral angle with
2203 # mirrored torque would also be fixed
2204 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2205 factor = 1.
2206 else:
2207 factor = 0.
2208 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \
2209 np.linalg.norm(p[1] - p0)
2210 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \
2211 np.linalg.norm(p[2] - p0)
2212 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2214 def index_shuffle(self, atoms, ind):
2215 # See docstring of superclass
2216 indices = []
2217 for new, old in slice2enlist(ind, len(atoms)):
2218 if old in self.indices:
2219 indices.append(new)
2220 if len(indices) == 0:
2221 raise IndexError('All indices in MirrorTorque not part of slice')
2222 self.indices = np.asarray(indices, int)
2224 def __repr__(self):
2225 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (
2226 self.indices[0], self.indices[1], self.indices[2],
2227 self.indices[3], self.max_angle, self.min_angle, self.fmax)
2229 def todict(self):
2230 return {'name': 'MirrorTorque',
2231 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2232 'a3': self.indices[2], 'a4': self.indices[3],
2233 'max_angle': self.max_angle,
2234 'min_angle': self.min_angle, 'fmax': self.fmax}}
2237class Filter(FilterOld):
2238 @deprecated('Import Filter from ase.filters')
2239 def __init__(self, *args, **kwargs):
2240 """
2241 .. deprecated:: 3.23.0
2242 Import ``Filter`` from :mod:`ase.filters`
2243 """
2244 super().__init__(*args, **kwargs)
2247class StrainFilter(StrainFilterOld):
2248 @deprecated('Import StrainFilter from ase.filters')
2249 def __init__(self, *args, **kwargs):
2250 """
2251 .. deprecated:: 3.23.0
2252 Import ``StrainFilter`` from :mod:`ase.filters`
2253 """
2254 super().__init__(*args, **kwargs)
2257class UnitCellFilter(UnitCellFilterOld):
2258 @deprecated('Import UnitCellFilter from ase.filters')
2259 def __init__(self, *args, **kwargs):
2260 """
2261 .. deprecated:: 3.23.0
2262 Import ``UnitCellFilter`` from :mod:`ase.filters`
2263 """
2264 super().__init__(*args, **kwargs)
2267class ExpCellFilter(ExpCellFilterOld):
2268 @deprecated('Import ExpCellFilter from ase.filters')
2269 def __init__(self, *args, **kwargs):
2270 """
2271 .. deprecated:: 3.23.0
2272 Import ``ExpCellFilter`` from :mod:`ase.filters`
2273 or use :class:`~ase.filters.FrechetCellFilter` for better
2274 convergence w.r.t. cell variables
2275 """
2276 super().__init__(*args, **kwargs)