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