Coverage for /builds/ase/ase/ase/constraints.py : 86.32%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1from warnings import warn
3import numpy as np
4from ase.calculators.calculator import PropertyNotImplementedError
5from ase.geometry import (find_mic, wrap_positions, get_distances_derivatives,
6 get_angles_derivatives, get_dihedrals_derivatives,
7 conditional_find_mic, get_angles, get_dihedrals)
8from ase.utils.parsemath import eval_expression
9from ase.stress import (full_3x3_to_voigt_6_stress,
10 voigt_6_to_full_3x3_stress)
12__all__ = [
13 'FixCartesian', 'FixBondLength', 'FixedMode',
14 'FixAtoms', 'UnitCellFilter', 'ExpCellFilter',
15 'FixScaled', 'StrainFilter', 'FixCom', 'FixedPlane', 'Filter',
16 'FixConstraint', 'FixedLine', 'FixBondLengths', 'FixLinearTriatomic',
17 'FixInternals', 'Hookean', 'ExternalForce', 'MirrorForce', 'MirrorTorque',
18 "FixScaledParametricRelations", "FixCartesianParametricRelations"]
21def dict2constraint(dct):
22 if dct['name'] not in __all__:
23 raise ValueError
24 return globals()[dct['name']](**dct['kwargs'])
27def slice2enlist(s, n):
28 """Convert a slice object into a list of (new, old) tuples."""
29 if isinstance(s, slice):
30 return enumerate(range(*s.indices(n)))
31 return enumerate(s)
34def constrained_indices(atoms, only_include=None):
35 """Returns a list of indices for the atoms that are constrained
36 by a constraint that is applied. By setting only_include to a
37 specific type of constraint you can make it only look for that
38 given constraint.
39 """
40 indices = []
41 for constraint in atoms.constraints:
42 if only_include is not None:
43 if not isinstance(constraint, only_include):
44 continue
45 indices.extend(np.array(constraint.get_indices()))
46 return np.array(np.unique(indices))
49class FixConstraint:
50 """Base class for classes that fix one or more atoms in some way."""
52 def index_shuffle(self, atoms, ind):
53 """Change the indices.
55 When the ordering of the atoms in the Atoms object changes,
56 this method can be called to shuffle the indices of the
57 constraints.
59 ind -- List or tuple of indices.
61 """
62 raise NotImplementedError
64 def repeat(self, m, n):
65 """ basic method to multiply by m, needs to know the length
66 of the underlying atoms object for the assignment of
67 multiplied constraints to work.
68 """
69 msg = ("Repeat is not compatible with your atoms' constraints."
70 ' Use atoms.set_constraint() before calling repeat to '
71 'remove your constraints.')
72 raise NotImplementedError(msg)
74 def adjust_momenta(self, atoms, momenta):
75 """Adjusts momenta in identical manner to forces."""
76 self.adjust_forces(atoms, momenta)
78 def copy(self):
79 return dict2constraint(self.todict().copy())
82class IndexedConstraint(FixConstraint):
83 def __init__(self, indices=None, mask=None):
84 """Constrain chosen atoms.
86 Parameters
87 ----------
88 indices : sequence of int
89 Indices for those atoms that should be constrained.
90 mask : sequence of bool
91 One boolean per atom indicating if the atom should be
92 constrained or not.
93 """
95 if mask is not None:
96 if indices is not None:
97 raise ValueError('Use only one of "indices" and "mask".')
98 indices = mask
99 indices = np.atleast_1d(indices)
100 if np.ndim(indices) > 1:
101 raise ValueError('indices has wrong amount of dimensions. '
102 f'Got {np.ndim(indices)}, expected ndim <= 1')
104 if indices.dtype == bool:
105 indices = np.arange(len(indices))[indices]
106 elif len(indices) == 0:
107 indices = np.empty(0, dtype=int)
108 elif not np.issubdtype(indices.dtype, np.integer):
109 raise ValueError('Indices must be integers or boolean mask, '
110 f'not dtype={indices.dtype}')
112 if len(set(indices)) < len(indices):
113 raise ValueError(
114 'The indices array contains duplicates. '
115 'Perhaps you want to specify a mask instead, but '
116 'forgot the mask= keyword.')
118 self.index = indices
120 def index_shuffle(self, atoms, ind):
121 # See docstring of superclass
122 index = []
124 # Resolve negative indices:
125 actual_indices = set(np.arange(len(atoms))[self.index])
127 for new, old in slice2enlist(ind, len(atoms)):
128 if old in actual_indices:
129 index.append(new)
130 if len(index) == 0:
131 raise IndexError('All indices in FixAtoms not part of slice')
132 self.index = np.asarray(index, int)
133 # XXX make immutable
135 def get_indices(self):
136 return self.index.copy()
138 def repeat(self, m, n):
139 i0 = 0
140 natoms = 0
141 if isinstance(m, int):
142 m = (m, m, m)
143 index_new = []
144 for m2 in range(m[2]):
145 for m1 in range(m[1]):
146 for m0 in range(m[0]):
147 i1 = i0 + n
148 index_new += [i + natoms for i in self.index]
149 i0 = i1
150 natoms += n
151 self.index = np.asarray(index_new, int)
152 # XXX make immutable
153 return self
155 def delete_atoms(self, indices, natoms):
156 """Removes atoms from the index array, if present.
158 Required for removing atoms with existing constraint.
159 """
161 i = np.zeros(natoms, int) - 1
162 new = np.delete(np.arange(natoms), indices)
163 i[new] = np.arange(len(new))
164 index = i[self.index]
165 self.index = index[index >= 0]
166 # XXX make immutable
167 if len(self.index) == 0:
168 return None
169 return self
172class FixAtoms(IndexedConstraint):
173 """Fix chosen atoms.
175 Examples
176 --------
177 Fix all Copper atoms:
179 >>> mask = (atoms.symbols == 'Cu')
180 >>> c = FixAtoms(mask=mask)
181 >>> atoms.set_constraint(c)
183 Fix all atoms with z-coordinate less than 1.0 Angstrom:
185 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)
186 >>> atoms.set_constraint(c)
187 """
189 def get_removed_dof(self, atoms):
190 return 3 * len(self.index)
192 def adjust_positions(self, atoms, new):
193 new[self.index] = atoms.positions[self.index]
195 def adjust_forces(self, atoms, forces):
196 forces[self.index] = 0.0
198 def __repr__(self):
199 clsname = type(self).__name__
200 indices = ints2string(self.index)
201 return f'{clsname}(indices={indices})'
203 def todict(self):
204 return {'name': 'FixAtoms',
205 'kwargs': {'indices': self.index.tolist()}}
208class FixCom(FixConstraint):
209 """Constraint class for fixing the center of mass.
211 References
213 https://pubs.acs.org/doi/abs/10.1021/jp9722824
215 """
217 def get_removed_dof(self, atoms):
218 return 3
220 def adjust_positions(self, atoms, new):
221 masses = atoms.get_masses()
222 old_cm = atoms.get_center_of_mass()
223 new_cm = np.dot(masses, new) / masses.sum()
224 d = old_cm - new_cm
225 new += d
227 def adjust_forces(self, atoms, forces):
228 m = atoms.get_masses()
229 mm = np.tile(m, (3, 1)).T
230 lb = np.sum(mm * forces, axis=0) / sum(m**2)
231 forces -= mm * lb
233 def todict(self):
234 return {'name': 'FixCom',
235 'kwargs': {}}
238def ints2string(x, threshold=None):
239 """Convert ndarray of ints to string."""
240 if threshold is None or len(x) <= threshold:
241 return str(x.tolist())
242 return str(x[:threshold].tolist())[:-1] + ', ...]'
245class FixBondLengths(FixConstraint):
246 maxiter = 500
248 def __init__(self, pairs, tolerance=1e-13,
249 bondlengths=None, iterations=None):
250 """iterations:
251 Ignored"""
252 self.pairs = np.asarray(pairs)
253 self.tolerance = tolerance
254 self.bondlengths = bondlengths
256 def get_removed_dof(self, atoms):
257 return len(self.pairs)
259 def adjust_positions(self, atoms, new):
260 old = atoms.positions
261 masses = atoms.get_masses()
263 if self.bondlengths is None:
264 self.bondlengths = self.initialize_bond_lengths(atoms)
266 for i in range(self.maxiter):
267 converged = True
268 for j, ab in enumerate(self.pairs):
269 a = ab[0]
270 b = ab[1]
271 cd = self.bondlengths[j]
272 r0 = old[a] - old[b]
273 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
274 d1 = new[a] - new[b] - r0 + d0
275 m = 1 / (1 / masses[a] + 1 / masses[b])
276 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)
277 if abs(x) > self.tolerance:
278 new[a] += x * m / masses[a] * d0
279 new[b] -= x * m / masses[b] * d0
280 converged = False
281 if converged:
282 break
283 else:
284 raise RuntimeError('Did not converge')
286 def adjust_momenta(self, atoms, p):
287 old = atoms.positions
288 masses = atoms.get_masses()
290 if self.bondlengths is None:
291 self.bondlengths = self.initialize_bond_lengths(atoms)
293 for i in range(self.maxiter):
294 converged = True
295 for j, ab in enumerate(self.pairs):
296 a = ab[0]
297 b = ab[1]
298 cd = self.bondlengths[j]
299 d = old[a] - old[b]
300 d, _ = find_mic(d, atoms.cell, atoms.pbc)
301 dv = p[a] / masses[a] - p[b] / masses[b]
302 m = 1 / (1 / masses[a] + 1 / masses[b])
303 x = -np.dot(dv, d) / cd**2
304 if abs(x) > self.tolerance:
305 p[a] += x * m * d
306 p[b] -= x * m * d
307 converged = False
308 if converged:
309 break
310 else:
311 raise RuntimeError('Did not converge')
313 def adjust_forces(self, atoms, forces):
314 self.constraint_forces = -forces
315 self.adjust_momenta(atoms, forces)
316 self.constraint_forces += forces
318 def initialize_bond_lengths(self, atoms):
319 bondlengths = np.zeros(len(self.pairs))
321 for i, ab in enumerate(self.pairs):
322 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)
324 return bondlengths
326 def get_indices(self):
327 return np.unique(self.pairs.ravel())
329 def todict(self):
330 return {'name': 'FixBondLengths',
331 'kwargs': {'pairs': self.pairs.tolist(),
332 'tolerance': self.tolerance}}
334 def index_shuffle(self, atoms, ind):
335 """Shuffle the indices of the two atoms in this constraint"""
336 map = np.zeros(len(atoms), int)
337 map[ind] = 1
338 n = map.sum()
339 map[:] = -1
340 map[ind] = range(n)
341 pairs = map[self.pairs]
342 self.pairs = pairs[(pairs != -1).all(1)]
343 if len(self.pairs) == 0:
344 raise IndexError('Constraint not part of slice')
347def FixBondLength(a1, a2):
348 """Fix distance between atoms with indices a1 and a2."""
349 return FixBondLengths([(a1, a2)])
352class FixLinearTriatomic(FixConstraint):
353 """Holonomic constraints for rigid linear triatomic molecules."""
355 def __init__(self, triples):
356 """Apply RATTLE-type bond constraints between outer atoms n and m
357 and linear vectorial constraints to the position of central
358 atoms o to fix the geometry of linear triatomic molecules of the
359 type:
361 n--o--m
363 Parameters:
365 triples: list
366 Indices of the atoms forming the linear molecules to constrain
367 as triples. Sequence should be (n, o, m) or (m, o, n).
369 When using these constraints in molecular dynamics or structure
370 optimizations, atomic forces need to be redistributed within a
371 triple. The function redistribute_forces_optimization implements
372 the redistribution of forces for structure optimization, while
373 the function redistribute_forces_md implements the redistribution
374 for molecular dynamics.
376 References:
378 Ciccotti et al. Molecular Physics 47 (1982)
379 :doi:`10.1080/00268978200100942`
380 """
381 self.triples = np.asarray(triples)
382 if self.triples.shape[1] != 3:
383 raise ValueError('"triples" has wrong size')
384 self.bondlengths = None
386 def get_removed_dof(self, atoms):
387 return 4 * len(self.triples)
389 @property
390 def n_ind(self):
391 return self.triples[:, 0]
393 @property
394 def m_ind(self):
395 return self.triples[:, 2]
397 @property
398 def o_ind(self):
399 return self.triples[:, 1]
401 def initialize(self, atoms):
402 masses = atoms.get_masses()
403 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)
405 self.bondlengths = self.initialize_bond_lengths(atoms)
406 self.bondlengths_nm = self.bondlengths.sum(axis=1)
408 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]
409 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +
410 C1[:, 1] ** 2 * self.mass_n * self.mass_o +
411 self.mass_n * self.mass_m)
412 C2 = C1 / C2[:, None]
413 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]
414 C3 = C2 * self.mass_o[:, None] * C3[:, None]
415 C3[:, 1] *= -1
416 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T
417 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)
418 C4 = C1 / C4[:, None]
420 self.C1 = C1
421 self.C2 = C2
422 self.C3 = C3
423 self.C4 = C4
425 def adjust_positions(self, atoms, new):
426 old = atoms.positions
427 new_n, new_m, new_o = self.get_slices(new)
429 if self.bondlengths is None:
430 self.initialize(atoms)
432 r0 = old[self.n_ind] - old[self.m_ind]
433 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
434 d1 = new_n - new_m - r0 + d0
435 a = np.einsum('ij,ij->i', d0, d0)
436 b = np.einsum('ij,ij->i', d1, d0)
437 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2
438 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))
439 g = g[:, None] * self.C3
440 new_n -= g[:, 0, None] * d0
441 new_m += g[:, 1, None] * d0
442 if np.allclose(d0, r0):
443 new_o = (self.C1[:, 0, None] * new_n
444 + self.C1[:, 1, None] * new_m)
445 else:
446 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)
447 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)
448 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2
449 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)
451 self.set_slices(new_n, new_m, new_o, new)
453 def adjust_momenta(self, atoms, p):
454 old = atoms.positions
455 p_n, p_m, p_o = self.get_slices(p)
457 if self.bondlengths is None:
458 self.initialize(atoms)
460 mass_nn = self.mass_n[:, None]
461 mass_mm = self.mass_m[:, None]
462 mass_oo = self.mass_o[:, None]
464 d = old[self.n_ind] - old[self.m_ind]
465 d, _ = find_mic(d, atoms.cell, atoms.pbc)
466 dv = p_n / mass_nn - p_m / mass_mm
467 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2
468 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]
469 p_n -= k[:, 0, None] * mass_nn * d
470 p_m += k[:, 1, None] * mass_mm * d
471 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +
472 self.C1[:, 1, None] * p_m / mass_mm))
474 self.set_slices(p_n, p_m, p_o, p)
476 def adjust_forces(self, atoms, forces):
478 if self.bondlengths is None:
479 self.initialize(atoms)
481 A = self.C4 * np.diff(self.C1)
482 A[:, 0] *= -1
483 A -= 1
484 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]
485 A /= (A.sum(axis=1))[:, None]
487 self.constraint_forces = -forces
488 old = atoms.positions
490 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)
492 d = old[self.n_ind] - old[self.m_ind]
493 d, _ = find_mic(d, atoms.cell, atoms.pbc)
494 df = fr_n - fr_m
495 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2
496 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]
497 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]
498 forces[self.o_ind] = fr_o + k[:, None] * d * B
500 self.constraint_forces += forces
502 def redistribute_forces_optimization(self, forces):
503 """Redistribute forces within a triple when performing structure
504 optimizations.
506 The redistributed forces needs to be further adjusted using the
507 appropriate Lagrange multipliers as implemented in adjust_forces."""
508 forces_n, forces_m, forces_o = self.get_slices(forces)
509 C1_1 = self.C1[:, 0, None]
510 C1_2 = self.C1[:, 1, None]
511 C4_1 = self.C4[:, 0, None]
512 C4_2 = self.C4[:, 1, None]
514 fr_n = ((1 - C4_1 * C1_1) * forces_n -
515 C4_1 * (C1_2 * forces_m - forces_o))
516 fr_m = ((1 - C4_2 * C1_2) * forces_m -
517 C4_2 * (C1_1 * forces_n - forces_o))
518 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +
519 C4_1 * forces_n + C4_2 * forces_m)
521 return fr_n, fr_m, fr_o
523 def redistribute_forces_md(self, atoms, forces, rand=False):
524 """Redistribute forces within a triple when performing molecular
525 dynamics.
527 When rand=True, use the equations for random force terms, as
528 used e.g. by Langevin dynamics, otherwise apply the standard
529 equations for deterministic forces (see Ciccotti et al. Molecular
530 Physics 47 (1982))."""
531 if self.bondlengths is None:
532 self.initialize(atoms)
533 forces_n, forces_m, forces_o = self.get_slices(forces)
534 C1_1 = self.C1[:, 0, None]
535 C1_2 = self.C1[:, 1, None]
536 C2_1 = self.C2[:, 0, None]
537 C2_2 = self.C2[:, 1, None]
538 mass_nn = self.mass_n[:, None]
539 mass_mm = self.mass_m[:, None]
540 mass_oo = self.mass_o[:, None]
541 if rand:
542 mr1 = (mass_mm / mass_nn) ** 0.5
543 mr2 = (mass_oo / mass_nn) ** 0.5
544 mr3 = (mass_nn / mass_mm) ** 0.5
545 mr4 = (mass_oo / mass_mm) ** 0.5
546 else:
547 mr1 = 1.0
548 mr2 = 1.0
549 mr3 = 1.0
550 mr4 = 1.0
552 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -
553 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -
554 mr2 * mass_mm * mass_nn * forces_o))
556 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -
557 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -
558 mr4 * mass_mm * mass_nn * forces_o))
560 self.set_slices(fr_n, fr_m, 0.0, forces)
562 def get_slices(self, a):
563 a_n = a[self.n_ind]
564 a_m = a[self.m_ind]
565 a_o = a[self.o_ind]
567 return a_n, a_m, a_o
569 def set_slices(self, a_n, a_m, a_o, a):
570 a[self.n_ind] = a_n
571 a[self.m_ind] = a_m
572 a[self.o_ind] = a_o
574 def initialize_bond_lengths(self, atoms):
575 bondlengths = np.zeros((len(self.triples), 2))
577 for i in range(len(self.triples)):
578 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],
579 self.o_ind[i], mic=True)
580 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],
581 self.m_ind[i], mic=True)
583 return bondlengths
585 def get_indices(self):
586 return np.unique(self.triples.ravel())
588 def todict(self):
589 return {'name': 'FixLinearTriatomic',
590 'kwargs': {'triples': self.triples.tolist()}}
592 def index_shuffle(self, atoms, ind):
593 """Shuffle the indices of the three atoms in this constraint"""
594 map = np.zeros(len(atoms), int)
595 map[ind] = 1
596 n = map.sum()
597 map[:] = -1
598 map[ind] = range(n)
599 triples = map[self.triples]
600 self.triples = triples[(triples != -1).all(1)]
601 if len(self.triples) == 0:
602 raise IndexError('Constraint not part of slice')
605class FixedMode(FixConstraint):
606 """Constrain atoms to move along directions orthogonal to
607 a given mode only. Initialize with a mode, such as one produced by
608 ase.vibrations.Vibrations.get_mode()."""
610 def __init__(self, mode):
611 mode = np.asarray(mode)
612 self.mode = (mode / np.sqrt((mode**2).sum())).reshape(-1)
614 def get_removed_dof(self, atoms):
615 return len(atoms)
617 def adjust_positions(self, atoms, newpositions):
618 newpositions = newpositions.ravel()
619 oldpositions = atoms.positions.ravel()
620 step = newpositions - oldpositions
621 newpositions -= self.mode * np.dot(step, self.mode)
623 def adjust_forces(self, atoms, forces):
624 forces = forces.ravel()
625 forces -= self.mode * np.dot(forces, self.mode)
627 def index_shuffle(self, atoms, ind):
628 eps = 1e-12
629 mode = self.mode.reshape(-1, 3)
630 excluded = np.ones(len(mode), dtype=bool)
631 excluded[ind] = False
632 if (abs(mode[excluded]) > eps).any():
633 raise IndexError('All nonzero parts of mode not in slice')
634 self.mode = mode[ind].ravel()
636 def get_indices(self):
637 # This function will never properly work because it works on all
638 # atoms and it has no idea how to tell how many atoms it is
639 # attached to. If it is being used, surely the user knows
640 # everything is being constrained.
641 return []
643 def todict(self):
644 return {'name': 'FixedMode',
645 'kwargs': {'mode': self.mode.tolist()}}
647 def __repr__(self):
648 return 'FixedMode(%s)' % self.mode.tolist()
651def _normalize(direction):
652 if np.shape(direction) != (3,):
653 raise ValueError("len(direction) is {len(direction)}. Has to be 3")
655 direction = np.asarray(direction) / np.linalg.norm(direction)
656 return direction
659class FixedPlane(IndexedConstraint):
660 """
661 Constraint object for fixing chosen atoms to only move in a plane.
663 The plane is defined by its normal vector *direction*
664 """
666 def __init__(self, indices, direction):
667 """Constrain chosen atoms.
669 Parameters
670 ----------
671 indices : int or list of int
672 Index or indices for atoms that should be constrained
673 direction : list of 3 int
674 Direction of the normal vector
676 Examples
677 --------
678 Fix all Copper atoms to only move in the yz-plane:
680 >>> from ase.constraints import FixedPlane
681 >>> c = FixedPlane(
682 >>> indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
683 >>> direction=[1, 0, 0],
684 >>> )
685 >>> atoms.set_constraint(c)
687 or constrain a single atom with the index 0 to move in the xy-plane:
689 >>> c = FixedPlane(indices=0, direction=[0, 0, 1])
690 >>> atoms.set_constraint(c)
691 """
692 super().__init__(indices=indices)
693 self.dir = _normalize(direction)
695 def adjust_positions(self, atoms, newpositions):
696 step = newpositions[self.index] - atoms.positions[self.index]
697 newpositions[self.index] -= _projection(step, self.dir)
699 def adjust_forces(self, atoms, forces):
700 forces[self.index] -= _projection(forces[self.index], self.dir)
702 def get_removed_dof(self, atoms):
703 return len(self.index)
705 def todict(self):
706 return {
707 'name': 'FixedPlane',
708 'kwargs': {'indices': self.index.tolist(),
709 'direction': self.dir.tolist()}
710 }
712 def __repr__(self):
713 return f'FixedPlane(indices={self.index}, {self.dir.tolist()})'
716def _projection(vectors, direction):
717 dotprods = vectors @ direction
718 projection = direction[None, :] * dotprods[:, None]
719 return projection
722class FixedLine(IndexedConstraint):
723 """
724 Constrain an atom index or a list of atom indices to move on a line only.
726 The line is defined by its vector *direction*
727 """
729 def __init__(self, indices, direction):
730 """Constrain chosen atoms.
732 Parameters
733 ----------
734 indices : int or list of int
735 Index or indices for atoms that should be constrained
736 direction : list of 3 int
737 Direction of the vector defining the line
739 Examples
740 --------
741 Fix all Copper atoms to only move in the x-direction:
743 >>> from ase.constraints import FixedLine
744 >>> c = FixedLine(
745 >>> indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],
746 >>> direction=[1, 0, 0],
747 >>> )
748 >>> atoms.set_constraint(c)
750 or constrain a single atom with the index 0 to move in the z-direction:
752 >>> c = FixedLine(indices=0, direction=[0, 0, 1])
753 >>> atoms.set_constraint(c)
754 """
755 super().__init__(indices)
756 self.dir = _normalize(direction)
758 def adjust_positions(self, atoms, newpositions):
759 step = newpositions[self.index] - atoms.positions[self.index]
760 projection = _projection(step, self.dir)
761 newpositions[self.index] = atoms.positions[self.index] + projection
763 def adjust_forces(self, atoms, forces):
764 forces[self.index] = _projection(forces[self.index], self.dir)
766 def get_removed_dof(self, atoms):
767 return 2 * len(self.index)
769 def __repr__(self):
770 return f'FixedLine(indices={self.index}, {self.dir.tolist()})'
772 def todict(self):
773 return {
774 'name': 'FixedLine',
775 'kwargs': {'indices': self.index.tolist(),
776 'direction': self.dir.tolist()}
777 }
780class FixCartesian(IndexedConstraint):
781 'Fix an atom index *a* in the directions of the cartesian coordinates.'
783 def __init__(self, a, mask=(1, 1, 1)):
784 super().__init__(indices=a)
785 self.mask = ~np.asarray(mask, bool)
787 def get_removed_dof(self, atoms):
788 return (3 - self.mask.sum()) * len(self.index)
790 def adjust_positions(self, atoms, new):
791 step = new[self.index] - atoms.positions[self.index]
792 step *= self.mask[None, :]
793 new[self.index] = atoms.positions[self.index] + step
795 def adjust_forces(self, atoms, forces):
796 forces[self.index] *= self.mask[None, :]
798 def __repr__(self):
799 return 'FixCartesian(indices={}, mask={})'.format(
800 self.index.tolist(), list(~self.mask))
802 def todict(self):
803 return {'name': 'FixCartesian',
804 'kwargs': {'a': self.index.tolist(),
805 'mask': (~self.mask).tolist()}}
808class FixScaled(IndexedConstraint):
809 'Fix an atom index *a* in the directions of the unit vectors.'
811 def __init__(self, a, mask=(1, 1, 1), cell=None):
812 # XXX The unused cell keyword is there for compatibility
813 # with old trajectory files.
814 super().__init__(a)
815 self.mask = np.array(mask, bool)
817 def get_removed_dof(self, atoms):
818 return self.mask.sum() * len(self.index)
820 def adjust_positions(self, atoms, new):
821 cell = atoms.cell
822 scaled_old = cell.scaled_positions(atoms.positions[self.index])
823 scaled_new = cell.scaled_positions(new[self.index])
824 scaled_new[:, self.mask] = scaled_old[:, self.mask]
825 new[self.index] = cell.cartesian_positions(scaled_new)
827 def adjust_forces(self, atoms, forces):
828 # Forces are covariant to the coordinate transformation,
829 # use the inverse transformations
830 cell = atoms.cell
831 scaled_forces = cell.cartesian_positions(forces[self.index])
832 scaled_forces *= -(self.mask - 1)
833 forces[self.index] = cell.scaled_positions(scaled_forces)
835 def todict(self):
836 return {'name': 'FixScaled',
837 'kwargs': {'a': self.index.tolist(),
838 'mask': self.mask.tolist()}}
840 def __repr__(self):
841 return 'FixScaled({}, {})'.format(self.index.tolist(), self.mask)
844# TODO: Better interface might be to use dictionaries in place of very
845# nested lists/tuples
846class FixInternals(FixConstraint):
847 """Constraint object for fixing multiple internal coordinates.
849 Allows fixing bonds, angles, dihedrals as well as linear combinations
850 of bonds (bondcombos).
852 Please provide angular units in degrees using `angles_deg` and
853 `dihedrals_deg`.
854 Fixing planar angles is not supported at the moment.
855 """
857 def __init__(self, bonds=None, angles=None, dihedrals=None,
858 angles_deg=None, dihedrals_deg=None,
859 bondcombos=None,
860 mic=False, epsilon=1.e-7):
861 """
862 A constrained internal coordinate is defined as a nested list:
863 '[value, [atom indices]]'. The constraint is initialized with a list of
864 constrained internal coordinates, i.e. '[[value, [atom indices]], ...]'.
865 If 'value' is None, the current value of the coordinate is constrained.
867 Parameters
868 ----------
869 bonds: nested python list, optional
870 List with targetvalue and atom indices defining the fixed bonds,
871 i.e. [[targetvalue, [index0, index1]], ...]
873 angles_deg: nested python list, optional
874 List with targetvalue and atom indices defining the fixedangles,
875 i.e. [[targetvalue, [index0, index1, index3]], ...]
877 dihedrals_deg: nested python list, optional
878 List with targetvalue and atom indices defining the fixed dihedrals,
879 i.e. [[targetvalue, [index0, index1, index3]], ...]
881 bondcombos: nested python list, optional
882 List with targetvalue, atom indices and linear coefficient defining
883 the fixed linear combination of bonds,
884 i.e. [[targetvalue, [[index0, index1, coefficient_for_bond],
885 [index1, index2, coefficient_for_bond]]], ...]
887 mic: bool, optional, default: False
888 Minimum image convention.
890 epsilon: float, optional, default: 1e-7
891 Convergence criterion.
892 """
893 warn_msg = 'Please specify {} in degrees using the {} argument.'
894 if angles:
895 warn(FutureWarning(warn_msg.format('angles', 'angle_deg')))
896 angles = np.asarray(angles)
897 angles[:, 0] = angles[:, 0] / np.pi * 180
898 angles = angles.tolist()
899 else:
900 angles = angles_deg
901 if dihedrals:
902 warn(FutureWarning(warn_msg.format('dihedrals', 'dihedrals_deg')))
903 dihedrals = np.asarray(dihedrals)
904 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180
905 dihedrals = dihedrals.tolist()
906 else:
907 dihedrals = dihedrals_deg
909 self.bonds = bonds or []
910 self.angles = angles or []
911 self.dihedrals = dihedrals or []
912 self.bondcombos = bondcombos or []
913 self.mic = mic
914 self.epsilon = epsilon
916 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)
917 + len(self.bondcombos))
919 # Initialize these at run-time:
920 self.constraints = []
921 self.initialized = False
923 def get_removed_dof(self, atoms):
924 return self.n
926 def initialize(self, atoms):
927 if self.initialized:
928 return
929 masses = np.repeat(atoms.get_masses(), 3)
930 cell = None
931 pbc = None
932 if self.mic:
933 cell = atoms.cell
934 pbc = atoms.pbc
935 self.constraints = []
936 for data, ConstrClass in [(self.bonds, self.FixBondLengthAlt),
937 (self.angles, self.FixAngle),
938 (self.dihedrals, self.FixDihedral),
939 (self.bondcombos, self.FixBondCombo)]:
940 for datum in data:
941 targetvalue = datum[0]
942 if targetvalue is None: # set to current value
943 targetvalue = ConstrClass.get_value(atoms, datum[1],
944 self.mic)
945 constr = ConstrClass(targetvalue, datum[1], masses, cell, pbc)
946 self.constraints.append(constr)
947 self.initialized = True
949 @staticmethod
950 def get_bondcombo(atoms, indices, mic=False):
951 """Convenience function to return the value of the bondcombo coordinate
952 (linear combination of bond lengths) for the given Atoms object 'atoms'.
953 Example: Get the current value of the linear combination of two bond
954 lengths defined as `bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]`."""
955 c = sum(df[2] * atoms.get_distance(*df[:2], mic=mic) for df in indices)
956 return c
958 def get_subconstraint(self, atoms, definition):
959 """Get pointer to a specific subconstraint.
960 Identification by its definition via indices (and coefficients)."""
961 self.initialize(atoms)
962 for subconstr in self.constraints:
963 if type(definition[0]) == list: # identify Combo constraint
964 defin = [d + [c] for d, c in zip(subconstr.indices,
965 subconstr.coefs)]
966 if defin == definition:
967 return subconstr
968 else: # identify primitive constraints by their indices
969 if subconstr.indices == [definition]:
970 return subconstr
971 raise ValueError('Given `definition` not found on Atoms object.')
973 def shuffle_definitions(self, shuffle_dic, internal_type):
974 dfns = [] # definitions
975 for dfn in internal_type: # e.g. for bond in self.bonds
976 append = True
977 new_dfn = [dfn[0], list(dfn[1])]
978 for old in dfn[1]:
979 if old in shuffle_dic:
980 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]
981 else:
982 append = False
983 break
984 if append:
985 dfns.append(new_dfn)
986 return dfns
988 def shuffle_combos(self, shuffle_dic, internal_type):
989 dfns = [] # definitions
990 for dfn in internal_type: # i.e. for bondcombo in self.bondcombos
991 append = True
992 all_indices = [idx[0:-1] for idx in dfn[1]]
993 new_dfn = [dfn[0], list(dfn[1])]
994 for i, indices in enumerate(all_indices):
995 for old in indices:
996 if old in shuffle_dic:
997 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]
998 else:
999 append = False
1000 break
1001 if not append:
1002 break
1003 if append:
1004 dfns.append(new_dfn)
1005 return dfns
1007 def index_shuffle(self, atoms, ind):
1008 # See docstring of superclass
1009 self.initialize(atoms)
1010 shuffle_dic = dict(slice2enlist(ind, len(atoms)))
1011 shuffle_dic = {old: new for new, old in shuffle_dic.items()}
1012 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)
1013 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)
1014 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)
1015 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)
1016 self.initialized = False
1017 self.initialize(atoms)
1018 if len(self.constraints) == 0:
1019 raise IndexError('Constraint not part of slice')
1021 def get_indices(self):
1022 cons = []
1023 for dfn in self.bonds + self.dihedrals + self.angles:
1024 cons.extend(dfn[1])
1025 for dfn in self.bondcombos:
1026 for partial_dfn in dfn[1]:
1027 cons.extend(partial_dfn[0:-1]) # last index is the coefficient
1028 return list(set(cons))
1030 def todict(self):
1031 return {'name': 'FixInternals',
1032 'kwargs': {'bonds': self.bonds,
1033 'angles_deg': self.angles,
1034 'dihedrals_deg': self.dihedrals,
1035 'bondcombos': self.bondcombos,
1036 'mic': self.mic,
1037 'epsilon': self.epsilon}}
1039 def adjust_positions(self, atoms, newpos):
1040 self.initialize(atoms)
1041 for constraint in self.constraints:
1042 constraint.setup_jacobian(atoms.positions)
1043 for j in range(50):
1044 maxerr = 0.0
1045 for constraint in self.constraints:
1046 constraint.adjust_positions(atoms.positions, newpos)
1047 maxerr = max(abs(constraint.sigma), maxerr)
1048 if maxerr < self.epsilon:
1049 return
1050 msg = 'FixInternals.adjust_positions did not converge.'
1051 if any([constr.targetvalue > 175. or constr.targetvalue < 5. for constr
1052 in self.constraints if isinstance(constr, self.FixAngle)]):
1053 msg += (' This may be caused by an almost planar angle.'
1054 ' Support for planar angles would require the'
1055 ' implementation of ghost, i.e. dummy, atoms.'
1056 ' See issue #868.')
1057 raise ValueError(msg)
1059 def adjust_forces(self, atoms, forces):
1060 """Project out translations and rotations and all other constraints"""
1061 self.initialize(atoms)
1062 positions = atoms.positions
1063 N = len(forces)
1064 list2_constraints = list(np.zeros((6, N, 3)))
1065 tx, ty, tz, rx, ry, rz = list2_constraints
1067 list_constraints = [r.ravel() for r in list2_constraints]
1069 tx[:, 0] = 1.0
1070 ty[:, 1] = 1.0
1071 tz[:, 2] = 1.0
1072 ff = forces.ravel()
1074 # Calculate the center of mass
1075 center = positions.sum(axis=0) / N
1077 rx[:, 1] = -(positions[:, 2] - center[2])
1078 rx[:, 2] = positions[:, 1] - center[1]
1079 ry[:, 0] = positions[:, 2] - center[2]
1080 ry[:, 2] = -(positions[:, 0] - center[0])
1081 rz[:, 0] = -(positions[:, 1] - center[1])
1082 rz[:, 1] = positions[:, 0] - center[0]
1084 # Normalizing transl., rotat. constraints
1085 for r in list2_constraints:
1086 r /= np.linalg.norm(r.ravel())
1088 # Add all angle, etc. constraint vectors
1089 for constraint in self.constraints:
1090 constraint.setup_jacobian(positions)
1091 constraint.adjust_forces(positions, forces)
1092 list_constraints.insert(0, constraint.jacobian)
1093 # QR DECOMPOSITION - GRAM SCHMIDT
1095 list_constraints = [r.ravel() for r in list_constraints]
1096 aa = np.column_stack(list_constraints)
1097 (aa, bb) = np.linalg.qr(aa)
1098 # Projection
1099 hh = []
1100 for i, constraint in enumerate(self.constraints):
1101 hh.append(aa[:, i] * np.row_stack(aa[:, i]))
1103 txx = aa[:, self.n] * np.row_stack(aa[:, self.n])
1104 tyy = aa[:, self.n + 1] * np.row_stack(aa[:, self.n + 1])
1105 tzz = aa[:, self.n + 2] * np.row_stack(aa[:, self.n + 2])
1106 rxx = aa[:, self.n + 3] * np.row_stack(aa[:, self.n + 3])
1107 ryy = aa[:, self.n + 4] * np.row_stack(aa[:, self.n + 4])
1108 rzz = aa[:, self.n + 5] * np.row_stack(aa[:, self.n + 5])
1109 T = txx + tyy + tzz + rxx + ryy + rzz
1110 for vec in hh:
1111 T += vec
1112 ff = np.dot(T, np.row_stack(ff))
1113 forces[:, :] -= np.dot(T, np.row_stack(ff)).reshape(-1, 3)
1115 def __repr__(self):
1116 constraints = [repr(constr) for constr in self.constraints]
1117 return f'FixInternals(_copy_init={constraints}, epsilon={self.epsilon})'
1119 # Classes for internal use in FixInternals
1120 class FixInternalsBase:
1121 """Base class for subclasses of FixInternals."""
1123 def __init__(self, targetvalue, indices, masses, cell, pbc):
1124 self.targetvalue = targetvalue # constant target value
1125 self.indices = [defin[0:-1] for defin in indices] # indices, defs
1126 self.coefs = np.asarray([defin[-1] for defin in indices])
1127 self.masses = masses
1128 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix
1129 self.sigma = 1. # difference between current and target value
1130 self.projected_force = None # helps optimizers scan along constr.
1131 self.cell = cell
1132 self.pbc = pbc
1134 def finalize_jacobian(self, pos, n_internals, n, derivs):
1135 """Populate jacobian with derivatives for `n_internals` defined
1136 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""
1137 jacobian = np.zeros((n_internals, *pos.shape))
1138 for i, idx in enumerate(self.indices):
1139 for j in range(n):
1140 jacobian[i, idx[j]] = derivs[i, j]
1141 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))
1142 return self.coefs @ jacobian
1144 def finalize_positions(self, newpos):
1145 jacobian = self.jacobian / self.masses
1146 lamda = -self.sigma / (jacobian @ self.get_jacobian(newpos))
1147 dnewpos = lamda * jacobian
1148 newpos += dnewpos.reshape(newpos.shape)
1150 def adjust_forces(self, positions, forces):
1151 self.projected_forces = ((self.jacobian @ forces.ravel())
1152 * self.jacobian)
1153 self.jacobian /= np.linalg.norm(self.jacobian)
1155 class FixBondCombo(FixInternalsBase):
1156 """Constraint subobject for fixing linear combination of bond lengths
1157 within FixInternals.
1159 sum_i( coef_i * bond_length_i ) = constant
1160 """
1162 def get_jacobian(self, pos):
1163 bondvectors = [pos[k] - pos[h] for h, k in self.indices]
1164 derivs = get_distances_derivatives(bondvectors, cell=self.cell,
1165 pbc=self.pbc)
1166 return self.finalize_jacobian(pos, len(bondvectors), 2, derivs)
1168 def setup_jacobian(self, pos):
1169 self.jacobian = self.get_jacobian(pos)
1171 def adjust_positions(self, oldpos, newpos):
1172 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]
1173 (_, ), (dists, ) = conditional_find_mic([bondvectors],
1174 cell=self.cell,
1175 pbc=self.pbc)
1176 value = self.coefs @ dists
1177 self.sigma = value - self.targetvalue
1178 self.finalize_positions(newpos)
1180 @staticmethod
1181 def get_value(atoms, indices, mic):
1182 return FixInternals.get_bondcombo(atoms, indices, mic)
1184 def __repr__(self):
1185 return (f'FixBondCombo({self.targetvalue}, {self.indices}, '
1186 '{self.coefs})')
1188 class FixBondLengthAlt(FixBondCombo):
1189 """Constraint subobject for fixing bond length within FixInternals.
1190 Fix distance between atoms with indices a1, a2."""
1192 def __init__(self, targetvalue, indices, masses, cell, pbc):
1193 if targetvalue <= 0.:
1194 raise ZeroDivisionError('Invalid targetvalue for fixed bond')
1195 indices = [list(indices) + [1.]] # bond definition with coef 1.
1196 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1198 @staticmethod
1199 def get_value(atoms, indices, mic):
1200 return atoms.get_distance(*indices, mic=mic)
1202 def __repr__(self):
1203 return f'FixBondLengthAlt({self.targetvalue}, {self.indices})'
1205 class FixAngle(FixInternalsBase):
1206 """Constraint subobject for fixing an angle within FixInternals.
1208 Convergence is potentially problematic for angles very close to
1209 0 or 180 degrees as there is a singularity in the Cartesian derivative.
1210 Fixing planar angles is therefore not supported at the moment.
1211 """
1213 def __init__(self, targetvalue, indices, masses, cell, pbc):
1214 """Fix atom movement to construct a constant angle."""
1215 if targetvalue <= 0. or targetvalue >= 180.:
1216 raise ZeroDivisionError('Invalid targetvalue for fixed angle')
1217 indices = [list(indices) + [1.]] # angle definition with coef 1.
1218 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1220 def gather_vectors(self, pos):
1221 v0 = [pos[h] - pos[k] for h, k, l in self.indices]
1222 v1 = [pos[l] - pos[k] for h, k, l in self.indices]
1223 return v0, v1
1225 def get_jacobian(self, pos):
1226 v0, v1 = self.gather_vectors(pos)
1227 derivs = get_angles_derivatives(v0, v1, cell=self.cell,
1228 pbc=self.pbc)
1229 return self.finalize_jacobian(pos, len(v0), 3, derivs)
1231 def setup_jacobian(self, pos):
1232 self.jacobian = self.get_jacobian(pos)
1234 def adjust_positions(self, oldpos, newpos):
1235 v0, v1 = self.gather_vectors(newpos)
1236 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)
1237 self.sigma = value - self.targetvalue
1238 self.finalize_positions(newpos)
1240 @staticmethod
1241 def get_value(atoms, indices, mic):
1242 return atoms.get_angle(*indices, mic=mic)
1244 def __repr__(self):
1245 return f'FixAngle({self.targetvalue}, {self.indices})'
1247 class FixDihedral(FixInternalsBase):
1248 """Constraint subobject for fixing a dihedral angle within FixInternals.
1250 A dihedral becomes undefined when at least one of the inner two angles
1251 becomes planar. Make sure to avoid this situation.
1252 """
1254 def __init__(self, targetvalue, indices, masses, cell, pbc):
1255 indices = [list(indices) + [1.]] # dihedral def. with coef 1.
1256 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1258 def gather_vectors(self, pos):
1259 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]
1260 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]
1261 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]
1262 return v0, v1, v2
1264 def get_jacobian(self, pos):
1265 v0, v1, v2 = self.gather_vectors(pos)
1266 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,
1267 pbc=self.pbc)
1268 return self.finalize_jacobian(pos, len(v0), 4, derivs)
1270 def setup_jacobian(self, pos):
1271 self.jacobian = self.get_jacobian(pos)
1273 def adjust_positions(self, oldpos, newpos):
1274 v0, v1, v2 = self.gather_vectors(newpos)
1275 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1276 # apply minimum dihedral difference 'convention': (diff <= 180)
1277 self.sigma = (value - self.targetvalue + 180) % 360 - 180
1278 self.finalize_positions(newpos)
1280 @staticmethod
1281 def get_value(atoms, indices, mic):
1282 return atoms.get_dihedral(*indices, mic=mic)
1284 def __repr__(self):
1285 return f'FixDihedral({self.targetvalue}, {self.indices})'
1288class FixParametricRelations(FixConstraint):
1290 def __init__(
1291 self,
1292 indices,
1293 Jacobian,
1294 const_shift,
1295 params=None,
1296 eps=1e-12,
1297 use_cell=False,
1298 ):
1299 """Constrains the degrees of freedom to act in a reduced parameter
1300 space defined by the Jacobian
1302 These constraints are based off the work in:
1303 https://arxiv.org/abs/1908.01610
1305 The constraints linearly maps the full 3N degrees of freedom,
1306 where N is number of active lattice vectors/atoms onto a
1307 reduced subset of M free parameters, where M <= 3*N. The
1308 Jacobian matrix and constant shift vector map the full set of
1309 degrees of freedom onto the reduced parameter space.
1311 Currently the constraint is set up to handle either atomic
1312 positions or lattice vectors at one time, but not both. To do
1313 both simply add a two constraints for each set. This is done
1314 to keep the mathematics behind the operations separate.
1316 It would be possible to extend these constraints to allow
1317 non-linear transformations if functionality to update the
1318 Jacobian at each position update was included. This would
1319 require passing an update function evaluate it every time
1320 adjust_positions is callled. This is currently NOT supported,
1321 and there are no plans to implement it in the future.
1323 Args:
1324 indices (list of int): indices of the constrained atoms
1325 (if not None or empty then cell_indices must be None or Empty)
1326 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))):
1327 The Jacobian describing
1328 the parameter space transformation
1329 const_shift (np.ndarray(shape=(3*len(indices)))):
1330 A vector describing the constant term
1331 in the transformation not accounted for in the Jacobian
1332 params (list of str):
1333 parameters used in the parametric representation
1334 if None a list is generated based on the shape of the Jacobian
1335 eps (float): a small number to compare the similarity of
1336 numbers and set the precision used
1337 to generate the constraint expressions
1338 use_cell (bool): if True then act on the cell object
1340 """
1341 self.indices = np.array(indices)
1342 self.Jacobian = np.array(Jacobian)
1343 self.const_shift = np.array(const_shift)
1345 assert self.const_shift.shape[0] == 3 * len(self.indices)
1346 assert self.Jacobian.shape[0] == 3 * len(self.indices)
1348 self.eps = eps
1349 self.use_cell = use_cell
1351 if params is None:
1352 params = []
1353 if self.Jacobian.shape[1] > 0:
1354 int_fmt_str = "{:0" + \
1355 str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"
1356 for param_ind in range(self.Jacobian.shape[1]):
1357 params.append("param_" + int_fmt_str.format(param_ind))
1358 else:
1359 assert len(params) == self.Jacobian.shape[-1]
1361 self.params = params
1363 self.Jacobian_inv = np.linalg.inv(
1364 self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T
1366 @classmethod
1367 def from_expressions(cls, indices, params, expressions,
1368 eps=1e-12, use_cell=False):
1369 """Converts the expressions into a Jacobian Matrix/const_shift
1370 vector and constructs a FixParametricRelations constraint
1372 The expressions must be a list like object of size 3*N and
1373 elements must be ordered as:
1374 [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],
1375 where i, j, and k are the first, second and third
1376 component of the atomic position/lattice
1377 vector. Currently only linear operations are allowed to be
1378 included in the expressions so
1379 only terms like:
1380 - const * param_0
1381 - sqrt[const] * param_1
1382 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M
1383 where const is any real number and param_0, param_1, ..., param_M are
1384 the parameters passed in
1385 params, are allowed.
1387 For example, fractional atomic position constraints for wurtzite are:
1388 params = ["z1", "z2"]
1389 expressions = [
1390 "1.0/3.0", "2.0/3.0", "z1",
1391 "2.0/3.0", "1.0/3.0", "0.5 + z1",
1392 "1.0/3.0", "2.0/3.0", "z2",
1393 "2.0/3.0", "1.0/3.0", "0.5 + z2",
1394 ]
1396 For diamond are:
1397 params = []
1398 expressions = [
1399 "0.0", "0.0", "0.0",
1400 "0.25", "0.25", "0.25",
1401 ],
1403 and for stannite are
1404 params=["x4", "z4"]
1405 expressions = [
1406 "0.0", "0.0", "0.0",
1407 "0.0", "0.5", "0.5",
1408 "0.75", "0.25", "0.5",
1409 "0.25", "0.75", "0.5",
1410 "x4 + z4", "x4 + z4", "2*x4",
1411 "x4 - z4", "x4 - z4", "-2*x4",
1412 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",
1413 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",
1414 ]
1416 Args:
1417 indices (list of int): indices of the constrained atoms
1418 (if not None or empty then cell_indices must be None or Empty)
1419 params (list of str): parameters used in the
1420 parametric representation
1421 expressions (list of str): expressions used to convert from the
1422 parametric to the real space representation
1423 eps (float): a small number to compare the similarity of
1424 numbers and set the precision used
1425 to generate the constraint expressions
1426 use_cell (bool): if True then act on the cell object
1428 Returns:
1429 cls(
1430 indices,
1431 Jacobian generated from expressions,
1432 const_shift generated from expressions,
1433 params,
1434 eps-12,
1435 use_cell,
1436 )
1437 """
1438 Jacobian = np.zeros((3 * len(indices), len(params)))
1439 const_shift = np.zeros(3 * len(indices))
1441 for expr_ind, expression in enumerate(expressions):
1442 expression = expression.strip()
1444 # Convert subtraction to addition
1445 expression = expression.replace("-", "+(-1.0)*")
1446 if "+" == expression[0]:
1447 expression = expression[1:]
1448 elif "(+" == expression[:2]:
1449 expression = "(" + expression[2:]
1451 # Explicitly add leading zeros so when replacing param_1 with 0.0
1452 # param_11 does not become 0.01
1453 int_fmt_str = "{:0" + \
1454 str(int(np.ceil(np.log10(len(params) + 1)))) + "d}"
1456 param_dct = dict()
1457 param_map = dict()
1459 # Construct a standardized param template for A/B filling
1460 for param_ind, param in enumerate(params):
1461 param_str = "param_" + int_fmt_str.format(param_ind)
1462 param_map[param] = param_str
1463 param_dct[param_str] = 0.0
1465 # Replace the parameters according to the map
1466 # Sort by string length (long to short) to prevent cases like x11
1467 # becoming f"{param_map["x1"]}1"
1468 for param in sorted(params, key=lambda s: -1.0 * len(s)):
1469 expression = expression.replace(param, param_map[param])
1471 # Partial linearity check
1472 for express_sec in expression.split("+"):
1473 in_sec = [param in express_sec for param in param_dct]
1474 n_params_in_sec = len(np.where(np.array(in_sec))[0])
1475 if n_params_in_sec > 1:
1476 raise ValueError(
1477 "FixParametricRelations expressions must be linear.")
1479 const_shift[expr_ind] = float(
1480 eval_expression(expression, param_dct))
1482 for param_ind in range(len(params)):
1483 param_str = "param_" + int_fmt_str.format(param_ind)
1484 if param_str not in expression:
1485 Jacobian[expr_ind, param_ind] = 0.0
1486 continue
1487 param_dct[param_str] = 1.0
1488 test_1 = float(eval_expression(expression, param_dct))
1489 test_1 -= const_shift[expr_ind]
1490 Jacobian[expr_ind, param_ind] = test_1
1492 param_dct[param_str] = 2.0
1493 test_2 = float(eval_expression(expression, param_dct))
1494 test_2 -= const_shift[expr_ind]
1495 if abs(test_2 / test_1 - 2.0) > eps:
1496 raise ValueError(
1497 "FixParametricRelations expressions must be linear.")
1498 param_dct[param_str] = 0.0
1500 args = [
1501 indices,
1502 Jacobian,
1503 const_shift,
1504 params,
1505 eps,
1506 use_cell,
1507 ]
1508 if cls is FixScaledParametricRelations:
1509 args = args[:-1]
1510 return cls(*args)
1512 @property
1513 def expressions(self):
1514 """Generate the expressions represented by the current self.Jacobian
1515 and self.const_shift objects"""
1516 expressions = []
1517 per = int(round(-1 * np.log10(self.eps)))
1518 fmt_str = "{:." + str(per + 1) + "g}"
1519 for index, shift_val in enumerate(self.const_shift):
1520 exp = ""
1521 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(
1522 shift_val) > self.eps:
1523 exp += fmt_str.format(shift_val)
1525 param_exp = ""
1526 for param_index, jacob_val in enumerate(self.Jacobian[index]):
1527 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)
1528 if abs_jacob_val < self.eps:
1529 continue
1531 param = self.params[param_index]
1532 if param_exp or exp:
1533 if jacob_val > -1.0 * self.eps:
1534 param_exp += " + "
1535 else:
1536 param_exp += " - "
1537 elif (not exp) and (not param_exp) and (
1538 jacob_val < -1.0 * self.eps):
1539 param_exp += "-"
1541 if np.abs(abs_jacob_val - 1.0) <= self.eps:
1542 param_exp += "{:s}".format(param)
1543 else:
1544 param_exp += (fmt_str +
1545 "*{:s}").format(abs_jacob_val, param)
1547 exp += param_exp
1549 expressions.append(exp)
1550 return np.array(expressions).reshape((-1, 3))
1552 def todict(self):
1553 """Create a dictionary representation of the constraint"""
1554 return {
1555 "name": type(self).__name__,
1556 "kwargs": {
1557 "indices": self.indices,
1558 "params": self.params,
1559 "Jacobian": self.Jacobian,
1560 "const_shift": self.const_shift,
1561 "eps": self.eps,
1562 "use_cell": self.use_cell,
1563 }
1564 }
1566 def __repr__(self):
1567 """The str representation of the constraint"""
1568 if len(self.indices) > 1:
1569 indices_str = "[{:d}, ..., {:d}]".format(
1570 self.indices[0], self.indices[-1])
1571 else:
1572 indices_str = "[{:d}]".format(self.indices[0])
1574 if len(self.params) > 1:
1575 params_str = "[{:s}, ..., {:s}]".format(
1576 self.params[0], self.params[-1])
1577 elif len(self.params) == 1:
1578 params_str = "[{:s}]".format(self.params[0])
1579 else:
1580 params_str = "[]"
1582 return '{:s}({:s}, {:s}, ..., {:e})'.format(
1583 type(self).__name__,
1584 indices_str,
1585 params_str,
1586 self.eps
1587 )
1590class FixScaledParametricRelations(FixParametricRelations):
1592 def __init__(
1593 self,
1594 indices,
1595 Jacobian,
1596 const_shift,
1597 params=None,
1598 eps=1e-12,
1599 ):
1600 """The fractional coordinate version of FixParametricRelations
1602 All arguments are the same, but since this is for fractional
1603 coordinates use_cell is false"""
1604 super(FixScaledParametricRelations, self).__init__(
1605 indices,
1606 Jacobian,
1607 const_shift,
1608 params,
1609 eps,
1610 False,
1611 )
1613 def adjust_contravariant(self, cell, vecs, B):
1614 """Adjust the values of a set of vectors that are contravariant
1615 with the unit transformation"""
1616 scaled = cell.scaled_positions(vecs).flatten()
1617 scaled = self.Jacobian_inv @ (scaled - B)
1618 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))
1620 return cell.cartesian_positions(scaled)
1622 def adjust_positions(self, atoms, positions):
1623 """Adjust positions of the atoms to match the constraints"""
1624 positions[self.indices] = self.adjust_contravariant(
1625 atoms.cell,
1626 positions[self.indices],
1627 self.const_shift,
1628 )
1629 positions[self.indices] = self.adjust_B(
1630 atoms.cell, positions[self.indices])
1632 def adjust_B(self, cell, positions):
1633 """Wraps the positions back to the unit cell and adjust B to
1634 keep track of this change"""
1635 fractional = cell.scaled_positions(positions)
1636 wrapped_fractional = (fractional % 1.0) % 1.0
1637 self.const_shift += np.round(wrapped_fractional - fractional).flatten()
1638 return cell.cartesian_positions(wrapped_fractional)
1640 def adjust_momenta(self, atoms, momenta):
1641 """Adjust momenta of the atoms to match the constraints"""
1642 momenta[self.indices] = self.adjust_contravariant(
1643 atoms.cell,
1644 momenta[self.indices],
1645 np.zeros(self.const_shift.shape),
1646 )
1648 def adjust_forces(self, atoms, forces):
1649 """Adjust forces of the atoms to match the constraints"""
1650 # Forces are coavarient to the coordinate transformation, use the
1651 # inverse transformations
1652 cart2frac_jacob = np.zeros(2 * (3 * len(atoms),))
1653 for i_atom in range(len(atoms)):
1654 cart2frac_jacob[3 * i_atom:3 * (i_atom + 1),
1655 3 * i_atom:3 * (i_atom + 1)] = atoms.cell.T
1657 jacobian = cart2frac_jacob @ self.Jacobian
1658 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T
1660 reduced_forces = jacobian.T @ forces.flatten()
1661 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)
1663 def todict(self):
1664 """Create a dictionary representation of the constraint"""
1665 dct = super(FixScaledParametricRelations, self).todict()
1666 del dct["kwargs"]["use_cell"]
1667 return dct
1670class FixCartesianParametricRelations(FixParametricRelations):
1672 def __init__(
1673 self,
1674 indices,
1675 Jacobian,
1676 const_shift,
1677 params=None,
1678 eps=1e-12,
1679 use_cell=False,
1680 ):
1681 """The Cartesian coordinate version of FixParametricRelations"""
1682 super(FixCartesianParametricRelations, self).__init__(
1683 indices,
1684 Jacobian,
1685 const_shift,
1686 params,
1687 eps,
1688 use_cell,
1689 )
1691 def adjust_contravariant(self, vecs, B):
1692 """Adjust the values of a set of vectors that are contravariant with
1693 the unit transformation"""
1694 vecs = self.Jacobian_inv @ (vecs.flatten() - B)
1695 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))
1696 return vecs
1698 def adjust_positions(self, atoms, positions):
1699 """Adjust positions of the atoms to match the constraints"""
1700 if self.use_cell:
1701 return
1702 positions[self.indices] = self.adjust_contravariant(
1703 positions[self.indices],
1704 self.const_shift,
1705 )
1707 def adjust_momenta(self, atoms, momenta):
1708 """Adjust momenta of the atoms to match the constraints"""
1709 if self.use_cell:
1710 return
1711 momenta[self.indices] = self.adjust_contravariant(
1712 momenta[self.indices],
1713 np.zeros(self.const_shift.shape),
1714 )
1716 def adjust_forces(self, atoms, forces):
1717 """Adjust forces of the atoms to match the constraints"""
1718 if self.use_cell:
1719 return
1721 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()
1722 forces[self.indices] = (self.Jacobian_inv.T @
1723 forces_reduced).reshape(-1, 3)
1725 def adjust_cell(self, atoms, cell):
1726 """Adjust the cell of the atoms to match the constraints"""
1727 if not self.use_cell:
1728 return
1729 cell[self.indices] = self.adjust_contravariant(
1730 cell[self.indices],
1731 np.zeros(self.const_shift.shape),
1732 )
1734 def adjust_stress(self, atoms, stress):
1735 """Adjust the stress of the atoms to match the constraints"""
1736 if not self.use_cell:
1737 return
1739 stress_3x3 = voigt_6_to_full_3x3_stress(stress)
1740 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()
1741 stress_3x3[self.indices] = (
1742 self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)
1744 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)
1747class Hookean(FixConstraint):
1748 """Applies a Hookean restorative force between a pair of atoms, an atom
1749 and a point, or an atom and a plane."""
1751 def __init__(self, a1, a2, k, rt=None):
1752 """Forces two atoms to stay close together by applying no force if
1753 they are below a threshold length, rt, and applying a Hookean
1754 restorative force when the distance between them exceeds rt. Can
1755 also be used to tether an atom to a fixed point in space or to a
1756 distance above a plane.
1758 a1 : int
1759 Index of atom 1
1760 a2 : one of three options
1761 1) index of atom 2
1762 2) a fixed point in cartesian space to which to tether a1
1763 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.
1764 k : float
1765 Hooke's law (spring) constant to apply when distance
1766 exceeds threshold_length. Units of eV A^-2.
1767 rt : float
1768 The threshold length below which there is no force. The
1769 length is 1) between two atoms, 2) between atom and point.
1770 This argument is not supplied in case 3. Units of A.
1772 If a plane is specified, the Hooke's law force is applied if the atom
1773 is on the normal side of the plane. For instance, the plane with
1774 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z
1775 intercept of +7 and a normal vector pointing in the +z direction.
1776 If the atom has z > 7, then a downward force would be applied of
1777 k * (atom.z - 7). The same plane with the normal vector pointing in
1778 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).
1780 References:
1782 Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014)
1783 https://link.springer.com/article/10.1007%2Fs11244-013-0161-8
1784 """
1786 if isinstance(a2, int):
1787 self._type = 'two atoms'
1788 self.indices = [a1, a2]
1789 elif len(a2) == 3:
1790 self._type = 'point'
1791 self.index = a1
1792 self.origin = np.array(a2)
1793 elif len(a2) == 4:
1794 self._type = 'plane'
1795 self.index = a1
1796 self.plane = a2
1797 else:
1798 raise RuntimeError('Unknown type for a2')
1799 self.threshold = rt
1800 self.spring = k
1802 def get_removed_dof(self, atoms):
1803 return 0
1805 def todict(self):
1806 dct = {'name': 'Hookean'}
1807 dct['kwargs'] = {'rt': self.threshold,
1808 'k': self.spring}
1809 if self._type == 'two atoms':
1810 dct['kwargs']['a1'] = self.indices[0]
1811 dct['kwargs']['a2'] = self.indices[1]
1812 elif self._type == 'point':
1813 dct['kwargs']['a1'] = self.index
1814 dct['kwargs']['a2'] = self.origin
1815 elif self._type == 'plane':
1816 dct['kwargs']['a1'] = self.index
1817 dct['kwargs']['a2'] = self.plane
1818 else:
1819 raise NotImplementedError('Bad type: %s' % self._type)
1820 return dct
1822 def adjust_positions(self, atoms, newpositions):
1823 pass
1825 def adjust_momenta(self, atoms, momenta):
1826 pass
1828 def adjust_forces(self, atoms, forces):
1829 positions = atoms.positions
1830 if self._type == 'plane':
1831 A, B, C, D = self.plane
1832 x, y, z = positions[self.index]
1833 d = ((A * x + B * y + C * z + D) /
1834 np.sqrt(A**2 + B**2 + C**2))
1835 if d < 0:
1836 return
1837 magnitude = self.spring * d
1838 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))
1839 forces[self.index] += direction * magnitude
1840 return
1841 if self._type == 'two atoms':
1842 p1, p2 = positions[self.indices]
1843 elif self._type == 'point':
1844 p1 = positions[self.index]
1845 p2 = self.origin
1846 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1847 bondlength = np.linalg.norm(displace)
1848 if bondlength > self.threshold:
1849 magnitude = self.spring * (bondlength - self.threshold)
1850 direction = displace / np.linalg.norm(displace)
1851 if self._type == 'two atoms':
1852 forces[self.indices[0]] += direction * magnitude
1853 forces[self.indices[1]] -= direction * magnitude
1854 else:
1855 forces[self.index] += direction * magnitude
1857 def adjust_potential_energy(self, atoms):
1858 """Returns the difference to the potential energy due to an active
1859 constraint. (That is, the quantity returned is to be added to the
1860 potential energy.)"""
1861 positions = atoms.positions
1862 if self._type == 'plane':
1863 A, B, C, D = self.plane
1864 x, y, z = positions[self.index]
1865 d = ((A * x + B * y + C * z + D) /
1866 np.sqrt(A**2 + B**2 + C**2))
1867 if d > 0:
1868 return 0.5 * self.spring * d**2
1869 else:
1870 return 0.
1871 if self._type == 'two atoms':
1872 p1, p2 = positions[self.indices]
1873 elif self._type == 'point':
1874 p1 = positions[self.index]
1875 p2 = self.origin
1876 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1877 bondlength = np.linalg.norm(displace)
1878 if bondlength > self.threshold:
1879 return 0.5 * self.spring * (bondlength - self.threshold)**2
1880 else:
1881 return 0.
1883 def get_indices(self):
1884 if self._type == 'two atoms':
1885 return self.indices
1886 elif self._type == 'point':
1887 return self.index
1888 elif self._type == 'plane':
1889 return self.index
1891 def index_shuffle(self, atoms, ind):
1892 # See docstring of superclass
1893 if self._type == 'two atoms':
1894 newa = [-1, -1] # Signal error
1895 for new, old in slice2enlist(ind, len(atoms)):
1896 for i, a in enumerate(self.indices):
1897 if old == a:
1898 newa[i] = new
1899 if newa[0] == -1 or newa[1] == -1:
1900 raise IndexError('Constraint not part of slice')
1901 self.indices = newa
1902 elif (self._type == 'point') or (self._type == 'plane'):
1903 newa = -1 # Signal error
1904 for new, old in slice2enlist(ind, len(atoms)):
1905 if old == self.index:
1906 newa = new
1907 break
1908 if newa == -1:
1909 raise IndexError('Constraint not part of slice')
1910 self.index = newa
1912 def __repr__(self):
1913 if self._type == 'two atoms':
1914 return 'Hookean(%d, %d)' % tuple(self.indices)
1915 elif self._type == 'point':
1916 return 'Hookean(%d) to cartesian' % self.index
1917 else:
1918 return 'Hookean(%d) to plane' % self.index
1921class ExternalForce(FixConstraint):
1922 """Constraint object for pulling two atoms apart by an external force.
1924 You can combine this constraint for example with FixBondLength but make
1925 sure that *ExternalForce* comes first in the list if there are overlaps
1926 between atom1-2 and atom3-4:
1928 >>> con1 = ExternalForce(atom1, atom2, f_ext)
1929 >>> con2 = FixBondLength(atom3, atom4)
1930 >>> atoms.set_constraint([con1, con2])
1932 see ase/test/external_force.py"""
1934 def __init__(self, a1, a2, f_ext):
1935 self.indices = [a1, a2]
1936 self.external_force = f_ext
1938 def get_removed_dof(self, atoms):
1939 return 0
1941 def adjust_positions(self, atoms, new):
1942 pass
1944 def adjust_forces(self, atoms, forces):
1945 dist = np.subtract.reduce(atoms.positions[self.indices])
1946 force = self.external_force * dist / np.linalg.norm(dist)
1947 forces[self.indices] += (force, -force)
1949 def adjust_potential_energy(self, atoms):
1950 dist = np.subtract.reduce(atoms.positions[self.indices])
1951 return -np.linalg.norm(dist) * self.external_force
1953 def index_shuffle(self, atoms, ind):
1954 """Shuffle the indices of the two atoms in this constraint"""
1955 newa = [-1, -1] # Signal error
1956 for new, old in slice2enlist(ind, len(atoms)):
1957 for i, a in enumerate(self.indices):
1958 if old == a:
1959 newa[i] = new
1960 if newa[0] == -1 or newa[1] == -1:
1961 raise IndexError('Constraint not part of slice')
1962 self.indices = newa
1964 def __repr__(self):
1965 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],
1966 self.indices[1],
1967 self.external_force)
1969 def todict(self):
1970 return {'name': 'ExternalForce',
1971 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
1972 'f_ext': self.external_force}}
1975class MirrorForce(FixConstraint):
1976 """Constraint object for mirroring the force between two atoms.
1978 This class is designed to find a transition state with the help of a
1979 single optimization. It can be used if the transition state belongs to a
1980 bond breaking reaction. First the given bond length will be fixed until
1981 all other degrees of freedom are optimized, then the forces of the two
1982 atoms will be mirrored to find the transition state. The mirror plane is
1983 perpendicular to the connecting line of the atoms. Transition states in
1984 dependence of the force can be obtained by stretching the molecule and
1985 fixing its total length with *FixBondLength* or by using *ExternalForce*
1986 during the optimization with *MirrorForce*.
1988 Parameters
1989 ----------
1990 a1: int
1991 First atom index.
1992 a2: int
1993 Second atom index.
1994 max_dist: float
1995 Upper limit of the bond length interval where the transition state
1996 can be found.
1997 min_dist: float
1998 Lower limit of the bond length interval where the transition state
1999 can be found.
2000 fmax: float
2001 Maximum force used for the optimization.
2003 Notes
2004 -----
2005 You can combine this constraint for example with FixBondLength but make
2006 sure that *MirrorForce* comes first in the list if there are overlaps
2007 between atom1-2 and atom3-4:
2009 >>> con1 = MirrorForce(atom1, atom2)
2010 >>> con2 = FixBondLength(atom3, atom4)
2011 >>> atoms.set_constraint([con1, con2])
2013 """
2015 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):
2016 self.indices = [a1, a2]
2017 self.min_dist = min_dist
2018 self.max_dist = max_dist
2019 self.fmax = fmax
2021 def adjust_positions(self, atoms, new):
2022 pass
2024 def adjust_forces(self, atoms, forces):
2025 dist = np.subtract.reduce(atoms.positions[self.indices])
2026 d = np.linalg.norm(dist)
2027 if (d < self.min_dist) or (d > self.max_dist):
2028 # Stop structure optimization
2029 forces[:] *= 0
2030 return
2031 dist /= d
2032 df = np.subtract.reduce(forces[self.indices])
2033 f = df.dot(dist)
2034 con_saved = atoms.constraints
2035 try:
2036 con = [con for con in con_saved
2037 if not isinstance(con, MirrorForce)]
2038 atoms.set_constraint(con)
2039 forces_copy = atoms.get_forces()
2040 finally:
2041 atoms.set_constraint(con_saved)
2042 df1 = -1 / 2. * f * dist
2043 forces_copy[self.indices] += (df1, -df1)
2044 # Check if forces would be converged if the bond with mirrored forces
2045 # would also be fixed
2046 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2047 factor = 1.
2048 else:
2049 factor = 0.
2050 df1 = -(1 + factor) / 2. * f * dist
2051 forces[self.indices] += (df1, -df1)
2053 def index_shuffle(self, atoms, ind):
2054 """Shuffle the indices of the two atoms in this constraint
2056 """
2057 newa = [-1, -1] # Signal error
2058 for new, old in slice2enlist(ind, len(atoms)):
2059 for i, a in enumerate(self.indices):
2060 if old == a:
2061 newa[i] = new
2062 if newa[0] == -1 or newa[1] == -1:
2063 raise IndexError('Constraint not part of slice')
2064 self.indices = newa
2066 def __repr__(self):
2067 return 'MirrorForce(%d, %d, %f, %f, %f)' % (
2068 self.indices[0], self.indices[1], self.max_dist, self.min_dist,
2069 self.fmax)
2071 def todict(self):
2072 return {'name': 'MirrorForce',
2073 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2074 'max_dist': self.max_dist,
2075 'min_dist': self.min_dist, 'fmax': self.fmax}}
2078class MirrorTorque(FixConstraint):
2079 """Constraint object for mirroring the torque acting on a dihedral
2080 angle defined by four atoms.
2082 This class is designed to find a transition state with the help of a
2083 single optimization. It can be used if the transition state belongs to a
2084 cis-trans-isomerization with a change of dihedral angle. First the given
2085 dihedral angle will be fixed until all other degrees of freedom are
2086 optimized, then the torque acting on the dihedral angle will be mirrored
2087 to find the transition state. Transition states in
2088 dependence of the force can be obtained by stretching the molecule and
2089 fixing its total length with *FixBondLength* or by using *ExternalForce*
2090 during the optimization with *MirrorTorque*.
2092 This constraint can be used to find
2093 transition states of cis-trans-isomerization.
2095 a1 a4
2096 | |
2097 a2 __ a3
2099 Parameters
2100 ----------
2101 a1: int
2102 First atom index.
2103 a2: int
2104 Second atom index.
2105 a3: int
2106 Third atom index.
2107 a4: int
2108 Fourth atom index.
2109 max_angle: float
2110 Upper limit of the dihedral angle interval where the transition state
2111 can be found.
2112 min_angle: float
2113 Lower limit of the dihedral angle interval where the transition state
2114 can be found.
2115 fmax: float
2116 Maximum force used for the optimization.
2118 Notes
2119 -----
2120 You can combine this constraint for example with FixBondLength but make
2121 sure that *MirrorTorque* comes first in the list if there are overlaps
2122 between atom1-4 and atom5-6:
2124 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)
2125 >>> con2 = FixBondLength(atom5, atom6)
2126 >>> atoms.set_constraint([con1, con2])
2128 """
2130 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,
2131 fmax=0.1):
2132 self.indices = [a1, a2, a3, a4]
2133 self.min_angle = min_angle
2134 self.max_angle = max_angle
2135 self.fmax = fmax
2137 def adjust_positions(self, atoms, new):
2138 pass
2140 def adjust_forces(self, atoms, forces):
2141 angle = atoms.get_dihedral(self.indices[0], self.indices[1],
2142 self.indices[2], self.indices[3])
2143 angle *= np.pi / 180.
2144 if (angle < self.min_angle) or (angle > self.max_angle):
2145 # Stop structure optimization
2146 forces[:] *= 0
2147 return
2148 p = atoms.positions[self.indices]
2149 f = forces[self.indices]
2151 f0 = (f[1] + f[2]) / 2.
2152 ff = f - f0
2153 p0 = (p[2] + p[1]) / 2.
2154 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)
2155 fff = ff - np.cross(m0, p - p0)
2156 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \
2157 (p[1] - p0).dot(p[1] - p0)
2158 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \
2159 (p[2] - p0).dot(p[2] - p0)
2160 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \
2161 np.linalg.norm(p[1] - p0)
2162 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \
2163 np.linalg.norm(p[2] - p0)
2164 omegap = omegap1 + omegap2
2165 con_saved = atoms.constraints
2166 try:
2167 con = [con for con in con_saved
2168 if not isinstance(con, MirrorTorque)]
2169 atoms.set_constraint(con)
2170 forces_copy = atoms.get_forces()
2171 finally:
2172 atoms.set_constraint(con_saved)
2173 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \
2174 np.linalg.norm(p[1] - p0)
2175 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \
2176 np.linalg.norm(p[2] - p0)
2177 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2178 # Check if forces would be converged if the dihedral angle with
2179 # mirrored torque would also be fixed
2180 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2181 factor = 1.
2182 else:
2183 factor = 0.
2184 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \
2185 np.linalg.norm(p[1] - p0)
2186 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \
2187 np.linalg.norm(p[2] - p0)
2188 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2190 def index_shuffle(self, atoms, ind):
2191 # See docstring of superclass
2192 indices = []
2193 for new, old in slice2enlist(ind, len(atoms)):
2194 if old in self.indices:
2195 indices.append(new)
2196 if len(indices) == 0:
2197 raise IndexError('All indices in MirrorTorque not part of slice')
2198 self.indices = np.asarray(indices, int)
2200 def __repr__(self):
2201 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (
2202 self.indices[0], self.indices[1], self.indices[2],
2203 self.indices[3], self.max_angle, self.min_angle, self.fmax)
2205 def todict(self):
2206 return {'name': 'MirrorTorque',
2207 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2208 'a3': self.indices[2], 'a4': self.indices[3],
2209 'max_angle': self.max_angle,
2210 'min_angle': self.min_angle, 'fmax': self.fmax}}
2213class Filter:
2214 """Subset filter class."""
2216 def __init__(self, atoms, indices=None, mask=None):
2217 """Filter atoms.
2219 This filter can be used to hide degrees of freedom in an Atoms
2220 object.
2222 Parameters
2223 ----------
2224 indices : list of int
2225 Indices for those atoms that should remain visible.
2226 mask : list of bool
2227 One boolean per atom indicating if the atom should remain
2228 visible or not.
2230 If a Trajectory tries to save this object, it will instead
2231 save the underlying Atoms object. To prevent this, override
2232 the iterimages method.
2233 """
2235 self.atoms = atoms
2236 self.constraints = []
2237 # Make self.info a reference to the underlying atoms' info dictionary.
2238 self.info = self.atoms.info
2240 if indices is None and mask is None:
2241 raise ValueError('Use "indices" or "mask".')
2242 if indices is not None and mask is not None:
2243 raise ValueError('Use only one of "indices" and "mask".')
2245 if mask is not None:
2246 self.index = np.asarray(mask, bool)
2247 self.n = self.index.sum()
2248 else:
2249 self.index = np.asarray(indices, int)
2250 self.n = len(self.index)
2252 def iterimages(self):
2253 # Present the real atoms object to Trajectory and friends
2254 return self.atoms.iterimages()
2256 def get_cell(self):
2257 """Returns the computational cell.
2259 The computational cell is the same as for the original system.
2260 """
2261 return self.atoms.get_cell()
2263 def get_pbc(self):
2264 """Returns the periodic boundary conditions.
2266 The boundary conditions are the same as for the original system.
2267 """
2268 return self.atoms.get_pbc()
2270 def get_positions(self):
2271 'Return the positions of the visible atoms.'
2272 return self.atoms.get_positions()[self.index]
2274 def set_positions(self, positions, **kwargs):
2275 'Set the positions of the visible atoms.'
2276 pos = self.atoms.get_positions()
2277 pos[self.index] = positions
2278 self.atoms.set_positions(pos, **kwargs)
2280 positions = property(get_positions, set_positions,
2281 doc='Positions of the atoms')
2283 def get_momenta(self):
2284 'Return the momenta of the visible atoms.'
2285 return self.atoms.get_momenta()[self.index]
2287 def set_momenta(self, momenta, **kwargs):
2288 'Set the momenta of the visible atoms.'
2289 mom = self.atoms.get_momenta()
2290 mom[self.index] = momenta
2291 self.atoms.set_momenta(mom, **kwargs)
2293 def get_atomic_numbers(self):
2294 'Return the atomic numbers of the visible atoms.'
2295 return self.atoms.get_atomic_numbers()[self.index]
2297 def set_atomic_numbers(self, atomic_numbers):
2298 'Set the atomic numbers of the visible atoms.'
2299 z = self.atoms.get_atomic_numbers()
2300 z[self.index] = atomic_numbers
2301 self.atoms.set_atomic_numbers(z)
2303 def get_tags(self):
2304 'Return the tags of the visible atoms.'
2305 return self.atoms.get_tags()[self.index]
2307 def set_tags(self, tags):
2308 'Set the tags of the visible atoms.'
2309 tg = self.atoms.get_tags()
2310 tg[self.index] = tags
2311 self.atoms.set_tags(tg)
2313 def get_forces(self, *args, **kwargs):
2314 return self.atoms.get_forces(*args, **kwargs)[self.index]
2316 def get_stress(self, *args, **kwargs):
2317 return self.atoms.get_stress(*args, **kwargs)
2319 def get_stresses(self, *args, **kwargs):
2320 return self.atoms.get_stresses(*args, **kwargs)[self.index]
2322 def get_masses(self):
2323 return self.atoms.get_masses()[self.index]
2325 def get_potential_energy(self, **kwargs):
2326 """Calculate potential energy.
2328 Returns the potential energy of the full system.
2329 """
2330 return self.atoms.get_potential_energy(**kwargs)
2332 def get_chemical_symbols(self):
2333 return self.atoms.get_chemical_symbols()
2335 def get_initial_magnetic_moments(self):
2336 return self.atoms.get_initial_magnetic_moments()
2338 def get_calculator(self):
2339 """Returns the calculator.
2341 WARNING: The calculator is unaware of this filter, and sees a
2342 different number of atoms.
2343 """
2344 return self.atoms.calc
2346 @property
2347 def calc(self):
2348 return self.atoms.calc
2350 def get_celldisp(self):
2351 return self.atoms.get_celldisp()
2353 def has(self, name):
2354 'Check for existence of array.'
2355 return self.atoms.has(name)
2357 def __len__(self):
2358 'Return the number of movable atoms.'
2359 return self.n
2361 def __getitem__(self, i):
2362 'Return an atom.'
2363 return self.atoms[self.index[i]]
2366class StrainFilter(Filter):
2367 """Modify the supercell while keeping the scaled positions fixed.
2369 Presents the strain of the supercell as the generalized positions,
2370 and the global stress tensor (times the volume) as the generalized
2371 force.
2373 This filter can be used to relax the unit cell until the stress is
2374 zero. If MDMin is used for this, the timestep (dt) to be used
2375 depends on the system size. 0.01/x where x is a typical dimension
2376 seems like a good choice.
2378 The stress and strain are presented as 6-vectors, the order of the
2379 components follow the standard engingeering practice: xx, yy, zz,
2380 yz, xz, xy.
2382 """
2384 def __init__(self, atoms, mask=None, include_ideal_gas=False):
2385 """Create a filter applying a homogeneous strain to a list of atoms.
2387 The first argument, atoms, is the atoms object.
2389 The optional second argument, mask, is a list of six booleans,
2390 indicating which of the six independent components of the
2391 strain that are allowed to become non-zero. It defaults to
2392 [1,1,1,1,1,1].
2394 """
2396 self.strain = np.zeros(6)
2397 self.include_ideal_gas = include_ideal_gas
2399 if mask is None:
2400 mask = np.ones(6)
2401 else:
2402 mask = np.array(mask)
2404 Filter.__init__(self, atoms, mask=mask)
2405 self.mask = mask
2406 self.origcell = atoms.get_cell()
2408 def get_positions(self):
2409 return self.strain.reshape((2, 3)).copy()
2411 def set_positions(self, new):
2412 new = new.ravel() * self.mask
2413 eps = np.array([[1.0 + new[0], 0.5 * new[5], 0.5 * new[4]],
2414 [0.5 * new[5], 1.0 + new[1], 0.5 * new[3]],
2415 [0.5 * new[4], 0.5 * new[3], 1.0 + new[2]]])
2417 self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True)
2418 self.strain[:] = new
2420 def get_forces(self, **kwargs):
2421 stress = self.atoms.get_stress(include_ideal_gas=self.include_ideal_gas)
2422 return -self.atoms.get_volume() * (stress * self.mask).reshape((2, 3))
2424 def has(self, x):
2425 return self.atoms.has(x)
2427 def __len__(self):
2428 return 2
2431class UnitCellFilter(Filter):
2432 """Modify the supercell and the atom positions. """
2434 def __init__(self, atoms, mask=None,
2435 cell_factor=None,
2436 hydrostatic_strain=False,
2437 constant_volume=False,
2438 scalar_pressure=0.0):
2439 """Create a filter that returns the atomic forces and unit cell
2440 stresses together, so they can simultaneously be minimized.
2442 The first argument, atoms, is the atoms object. The optional second
2443 argument, mask, is a list of booleans, indicating which of the six
2444 independent components of the strain are relaxed.
2446 - True = relax to zero
2447 - False = fixed, ignore this component
2449 Degrees of freedom are the positions in the original undeformed cell,
2450 plus the deformation tensor (extra 3 "atoms"). This gives forces
2451 consistent with numerical derivatives of the potential energy
2452 with respect to the cell degreees of freedom.
2454 For full details see:
2455 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
2456 Phys. Rev. B 59, 235 (1999)
2458 You can still use constraints on the atoms, e.g. FixAtoms, to control
2459 the relaxation of the atoms.
2461 >>> # this should be equivalent to the StrainFilter
2462 >>> atoms = Atoms(...)
2463 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
2464 >>> ucf = UnitCellFilter(atoms)
2466 You should not attach this UnitCellFilter object to a
2467 trajectory. Instead, create a trajectory for the atoms, and
2468 attach it to an optimizer like this:
2470 >>> atoms = Atoms(...)
2471 >>> ucf = UnitCellFilter(atoms)
2472 >>> qn = QuasiNewton(ucf)
2473 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
2474 >>> qn.attach(traj)
2475 >>> qn.run(fmax=0.05)
2477 Helpful conversion table:
2479 - 0.05 eV/A^3 = 8 GPA
2480 - 0.003 eV/A^3 = 0.48 GPa
2481 - 0.0006 eV/A^3 = 0.096 GPa
2482 - 0.0003 eV/A^3 = 0.048 GPa
2483 - 0.0001 eV/A^3 = 0.02 GPa
2485 Additional optional arguments:
2487 cell_factor: float (default float(len(atoms)))
2488 Factor by which deformation gradient is multiplied to put
2489 it on the same scale as the positions when assembling
2490 the combined position/cell vector. The stress contribution to
2491 the forces is scaled down by the same factor. This can be thought
2492 of as a very simple preconditioners. Default is number of atoms
2493 which gives approximately the correct scaling.
2495 hydrostatic_strain: bool (default False)
2496 Constrain the cell by only allowing hydrostatic deformation.
2497 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
2499 constant_volume: bool (default False)
2500 Project out the diagonal elements of the virial tensor to allow
2501 relaxations at constant volume, e.g. for mapping out an
2502 energy-volume curve. Note: this only approximately conserves
2503 the volume and breaks energy/force consistency so can only be
2504 used with optimizers that do require do a line minimisation
2505 (e.g. FIRE).
2507 scalar_pressure: float (default 0.0)
2508 Applied pressure to use for enthalpy pV term. As above, this
2509 breaks energy/force consistency.
2510 """
2512 Filter.__init__(self, atoms, indices=range(len(atoms)))
2513 self.atoms = atoms
2514 self.orig_cell = atoms.get_cell()
2515 self.stress = None
2517 if mask is None:
2518 mask = np.ones(6)
2519 mask = np.asarray(mask)
2520 if mask.shape == (6,):
2521 self.mask = voigt_6_to_full_3x3_stress(mask)
2522 elif mask.shape == (3, 3):
2523 self.mask = mask
2524 else:
2525 raise ValueError('shape of mask should be (3,3) or (6,)')
2527 if cell_factor is None:
2528 cell_factor = float(len(atoms))
2529 self.hydrostatic_strain = hydrostatic_strain
2530 self.constant_volume = constant_volume
2531 self.scalar_pressure = scalar_pressure
2532 self.cell_factor = cell_factor
2533 self.copy = self.atoms.copy
2534 self.arrays = self.atoms.arrays
2536 def deform_grad(self):
2537 return np.linalg.solve(self.orig_cell, self.atoms.cell).T
2539 def get_positions(self):
2540 """
2541 this returns an array with shape (natoms + 3,3).
2543 the first natoms rows are the positions of the atoms, the last
2544 three rows are the deformation tensor associated with the unit cell,
2545 scaled by self.cell_factor.
2546 """
2548 cur_deform_grad = self.deform_grad()
2549 natoms = len(self.atoms)
2550 pos = np.zeros((natoms + 3, 3))
2551 # UnitCellFilter's positions are the self.atoms.positions but without
2552 # the applied deformation gradient
2553 pos[:natoms] = np.linalg.solve(cur_deform_grad,
2554 self.atoms.positions.T).T
2555 # UnitCellFilter's cell DOFs are the deformation gradient times a
2556 # scaling factor
2557 pos[natoms:] = self.cell_factor * cur_deform_grad
2558 return pos
2560 def set_positions(self, new, **kwargs):
2561 """
2562 new is an array with shape (natoms+3,3).
2564 the first natoms rows are the positions of the atoms, the last
2565 three rows are the deformation tensor used to change the cell shape.
2567 the new cell is first set from original cell transformed by the new
2568 deformation gradient, then the positions are set with respect to the
2569 current cell by transforming them with the same deformation gradient
2570 """
2572 natoms = len(self.atoms)
2573 new_atom_positions = new[:natoms]
2574 new_deform_grad = new[natoms:] / self.cell_factor
2575 # Set the new cell from the original cell and the new
2576 # deformation gradient. Both current and final structures should
2577 # preserve symmetry, so if set_cell() calls FixSymmetry.adjust_cell(),
2578 # it should be OK
2579 self.atoms.set_cell(self.orig_cell @ new_deform_grad.T,
2580 scale_atoms=True)
2581 # Set the positions from the ones passed in (which are without the
2582 # deformation gradient applied) and the new deformation gradient.
2583 # This should also preserve symmetry, so if set_positions() calls
2584 # FixSymmetyr.adjust_positions(), it should be OK
2585 self.atoms.set_positions(new_atom_positions @ new_deform_grad.T,
2586 **kwargs)
2588 def get_potential_energy(self, force_consistent=True):
2589 """
2590 returns potential energy including enthalpy PV term.
2591 """
2592 atoms_energy = self.atoms.get_potential_energy(
2593 force_consistent=force_consistent)
2594 return atoms_energy + self.scalar_pressure * self.atoms.get_volume()
2596 def get_forces(self, **kwargs):
2597 """
2598 returns an array with shape (natoms+3,3) of the atomic forces
2599 and unit cell stresses.
2601 the first natoms rows are the forces on the atoms, the last
2602 three rows are the forces on the unit cell, which are
2603 computed from the stress tensor.
2604 """
2606 stress = self.atoms.get_stress(**kwargs)
2607 atoms_forces = self.atoms.get_forces(**kwargs)
2609 volume = self.atoms.get_volume()
2610 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
2611 np.diag([self.scalar_pressure] * 3))
2612 cur_deform_grad = self.deform_grad()
2613 atoms_forces = atoms_forces @ cur_deform_grad
2614 virial = np.linalg.solve(cur_deform_grad, virial.T).T
2616 if self.hydrostatic_strain:
2617 vtr = virial.trace()
2618 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
2620 # Zero out components corresponding to fixed lattice elements
2621 if (self.mask != 1.0).any():
2622 virial *= self.mask
2624 if self.constant_volume:
2625 vtr = virial.trace()
2626 np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0)
2628 natoms = len(self.atoms)
2629 forces = np.zeros((natoms + 3, 3))
2630 forces[:natoms] = atoms_forces
2631 forces[natoms:] = virial / self.cell_factor
2633 self.stress = -full_3x3_to_voigt_6_stress(virial) / volume
2634 return forces
2636 def get_stress(self):
2637 raise PropertyNotImplementedError
2639 def has(self, x):
2640 return self.atoms.has(x)
2642 def __len__(self):
2643 return (len(self.atoms) + 3)
2646class ExpCellFilter(UnitCellFilter):
2647 """Modify the supercell and the atom positions."""
2649 def __init__(self, atoms, mask=None,
2650 cell_factor=None,
2651 hydrostatic_strain=False,
2652 constant_volume=False,
2653 scalar_pressure=0.0):
2654 r"""Create a filter that returns the atomic forces and unit cell
2655 stresses together, so they can simultaneously be minimized.
2657 The first argument, atoms, is the atoms object. The optional second
2658 argument, mask, is a list of booleans, indicating which of the six
2659 independent components of the strain are relaxed.
2661 - True = relax to zero
2662 - False = fixed, ignore this component
2664 Degrees of freedom are the positions in the original undeformed cell,
2665 plus the log of the deformation tensor (extra 3 "atoms"). This gives
2666 forces consistent with numerical derivatives of the potential energy
2667 with respect to the cell degrees of freedom.
2669 For full details see:
2670 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
2671 Phys. Rev. B 59, 235 (1999)
2673 You can still use constraints on the atoms, e.g. FixAtoms, to control
2674 the relaxation of the atoms.
2676 >>> # this should be equivalent to the StrainFilter
2677 >>> atoms = Atoms(...)
2678 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
2679 >>> ecf = ExpCellFilter(atoms)
2681 You should not attach this ExpCellFilter object to a
2682 trajectory. Instead, create a trajectory for the atoms, and
2683 attach it to an optimizer like this:
2685 >>> atoms = Atoms(...)
2686 >>> ecf = ExpCellFilter(atoms)
2687 >>> qn = QuasiNewton(ecf)
2688 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
2689 >>> qn.attach(traj)
2690 >>> qn.run(fmax=0.05)
2692 Helpful conversion table:
2694 - 0.05 eV/A^3 = 8 GPA
2695 - 0.003 eV/A^3 = 0.48 GPa
2696 - 0.0006 eV/A^3 = 0.096 GPa
2697 - 0.0003 eV/A^3 = 0.048 GPa
2698 - 0.0001 eV/A^3 = 0.02 GPa
2700 Additional optional arguments:
2702 cell_factor: (DEPRECATED)
2703 Retained for backwards compatibility, but no longer used.
2705 hydrostatic_strain: bool (default False)
2706 Constrain the cell by only allowing hydrostatic deformation.
2707 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
2709 constant_volume: bool (default False)
2710 Project out the diagonal elements of the virial tensor to allow
2711 relaxations at constant volume, e.g. for mapping out an
2712 energy-volume curve.
2714 scalar_pressure: float (default 0.0)
2715 Applied pressure to use for enthalpy pV term. As above, this
2716 breaks energy/force consistency.
2718 Implementation details:
2720 The implementation is based on that of Christoph Ortner in JuLIP.jl:
2721 https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244
2723 We decompose the deformation gradient as
2725 F = exp(U) F0
2726 x = F * F0^{-1} z = exp(U) z
2728 If we write the energy as a function of U we can transform the
2729 stress associated with a perturbation V into a derivative using a
2730 linear map V -> L(U, V).
2732 \phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) .
2733 ( L(U, V) exp(-U) exp(U) z )
2735 where
2737 \nabla E(U) : V = [S exp(-U)'] : L(U,V)
2738 = L'(U, S exp(-U)') : V
2739 = L(U', S exp(-U)') : V
2740 = L(U, S exp(-U)) : V (provided U = U')
2742 where the : operator represents double contraction,
2743 i.e. A:B = trace(A'B), and
2745 F = deformation tensor - 3x3 matrix
2746 F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here
2747 U = cell degrees of freedom used here - 3x3 matrix
2748 V = perturbation to cell DoFs - 3x3 matrix
2749 v = perturbation to position DoFs
2750 x = atomic positions in deformed cell
2751 z = atomic positions in original cell
2752 \phi = potential energy
2753 S = stress tensor [3x3 matrix]
2754 L(U, V) = directional derivative of exp at U in direction V, i.e
2755 d/dt exp(U + t V)|_{t=0} = L(U, V)
2757 This means we can write
2759 d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V
2761 and therefore the contribution to the gradient of the energy is
2763 \nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij
2764 """
2766 Filter.__init__(self, atoms, indices=range(len(atoms)))
2767 UnitCellFilter.__init__(self, atoms, mask, cell_factor,
2768 hydrostatic_strain,
2769 constant_volume, scalar_pressure)
2770 if cell_factor is not None:
2771 warn("cell_factor is deprecated")
2772 self.cell_factor = 1.0
2774 # We defer the scipy import to avoid high immediate import overhead
2775 from scipy.linalg import expm, logm
2776 self.expm = expm
2777 self.logm = logm
2779 def get_positions(self):
2780 pos = UnitCellFilter.get_positions(self)
2781 natoms = len(self.atoms)
2782 pos[natoms:] = self.logm(self.deform_grad())
2783 return pos
2785 def set_positions(self, new, **kwargs):
2786 natoms = len(self.atoms)
2787 new2 = new.copy()
2788 new2[natoms:] = self.expm(new[natoms:])
2789 UnitCellFilter.set_positions(self, new2, **kwargs)
2791 def get_forces(self, **kwargs):
2792 forces = UnitCellFilter.get_forces(self, **kwargs)
2794 # forces on atoms are same as UnitCellFilter, we just
2795 # need to modify the stress contribution
2796 stress = self.atoms.get_stress(**kwargs)
2797 volume = self.atoms.get_volume()
2798 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
2799 np.diag([self.scalar_pressure] * 3))
2801 cur_deform_grad = self.deform_grad()
2802 cur_deform_grad_log = self.logm(cur_deform_grad)
2804 if self.hydrostatic_strain:
2805 vtr = virial.trace()
2806 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
2808 # Zero out components corresponding to fixed lattice elements
2809 if (self.mask != 1.0).any():
2810 virial *= self.mask
2812 deform_grad_log_force_naive = virial.copy()
2813 Y = np.zeros((6, 6))
2814 Y[0:3, 0:3] = cur_deform_grad_log
2815 Y[3:6, 3:6] = cur_deform_grad_log
2816 Y[0:3, 3:6] = - virial @ self.expm(-cur_deform_grad_log)
2817 deform_grad_log_force = -self.expm(Y)[0:3, 3:6]
2818 for (i1, i2) in [(0, 1), (0, 2), (1, 2)]:
2819 ff = 0.5 * (deform_grad_log_force[i1, i2] +
2820 deform_grad_log_force[i2, i1])
2821 deform_grad_log_force[i1, i2] = ff
2822 deform_grad_log_force[i2, i1] = ff
2824 # check for reasonable alignment between naive and
2825 # exact search directions
2826 all_are_equal = np.all(np.isclose(deform_grad_log_force,
2827 deform_grad_log_force_naive))
2828 if all_are_equal or \
2829 (np.sum(deform_grad_log_force * deform_grad_log_force_naive) /
2830 np.sqrt(np.sum(deform_grad_log_force**2) *
2831 np.sum(deform_grad_log_force_naive**2)) > 0.8):
2832 deform_grad_log_force = deform_grad_log_force_naive
2834 # Cauchy stress used for convergence testing
2835 convergence_crit_stress = -(virial / volume)
2836 if self.constant_volume:
2837 # apply constraint to force
2838 dglf_trace = deform_grad_log_force.trace()
2839 np.fill_diagonal(deform_grad_log_force,
2840 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
2841 # apply constraint to Cauchy stress used for convergence testing
2842 ccs_trace = convergence_crit_stress.trace()
2843 np.fill_diagonal(convergence_crit_stress,
2844 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
2846 # pack gradients into vector
2847 natoms = len(self.atoms)
2848 forces[natoms:] = deform_grad_log_force
2849 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
2850 return forces