Coverage for /builds/debichem-team/python-ase/ase/constraints.py : 83.50%
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 math import sqrt
2from warnings import warn
4import numpy as np
5from scipy.linalg import expm, logm
6from ase.calculators.calculator import PropertyNotImplementedError
7from ase.geometry import (find_mic, wrap_positions, get_distances_derivatives,
8 get_angles_derivatives, get_dihedrals_derivatives,
9 conditional_find_mic, get_angles, get_dihedrals)
10from ase.utils.parsemath import eval_expression
11from ase.stress import (full_3x3_to_voigt_6_stress,
12 voigt_6_to_full_3x3_stress)
14__all__ = [
15 'FixCartesian', 'FixBondLength', 'FixedMode',
16 'FixConstraintSingle', 'FixAtoms', 'UnitCellFilter', 'ExpCellFilter',
17 'FixScaled', 'StrainFilter', 'FixCom', 'FixedPlane', 'Filter',
18 'FixConstraint', 'FixedLine', 'FixBondLengths', 'FixLinearTriatomic',
19 'FixInternals', 'Hookean', 'ExternalForce', 'MirrorForce', 'MirrorTorque',
20 "FixScaledParametricRelations", "FixCartesianParametricRelations"]
23def dict2constraint(dct):
24 if dct['name'] not in __all__:
25 raise ValueError
26 return globals()[dct['name']](**dct['kwargs'])
29def slice2enlist(s, n):
30 """Convert a slice object into a list of (new, old) tuples."""
31 if isinstance(s, slice):
32 return enumerate(range(*s.indices(n)))
33 return enumerate(s)
36def constrained_indices(atoms, only_include=None):
37 """Returns a list of indices for the atoms that are constrained
38 by a constraint that is applied. By setting only_include to a
39 specific type of constraint you can make it only look for that
40 given constraint.
41 """
42 indices = []
43 for constraint in atoms.constraints:
44 if only_include is not None:
45 if not isinstance(constraint, only_include):
46 continue
47 indices.extend(np.array(constraint.get_indices()))
48 return np.array(np.unique(indices))
51class FixConstraint:
52 """Base class for classes that fix one or more atoms in some way."""
54 def index_shuffle(self, atoms, ind):
55 """Change the indices.
57 When the ordering of the atoms in the Atoms object changes,
58 this method can be called to shuffle the indices of the
59 constraints.
61 ind -- List or tuple of indices.
63 """
64 raise NotImplementedError
66 def repeat(self, m, n):
67 """ basic method to multiply by m, needs to know the length
68 of the underlying atoms object for the assignment of
69 multiplied constraints to work.
70 """
71 msg = ("Repeat is not compatible with your atoms' constraints."
72 ' Use atoms.set_constraint() before calling repeat to '
73 'remove your constraints.')
74 raise NotImplementedError(msg)
76 def adjust_momenta(self, atoms, momenta):
77 """Adjusts momenta in identical manner to forces."""
78 self.adjust_forces(atoms, momenta)
80 def copy(self):
81 return dict2constraint(self.todict().copy())
84class FixConstraintSingle(FixConstraint):
85 """Base class for classes that fix a single atom."""
87 def __init__(self, a):
88 self.a = a
90 def index_shuffle(self, atoms, ind):
91 """The atom index must be stored as self.a."""
92 newa = None # Signal error
93 if self.a < 0:
94 self.a += len(atoms)
95 for new, old in slice2enlist(ind, len(atoms)):
96 if old == self.a:
97 newa = new
98 break
99 if newa is None:
100 raise IndexError('Constraint not part of slice')
101 self.a = newa
103 def get_indices(self):
104 return [self.a]
107class FixAtoms(FixConstraint):
108 """Constraint object for fixing some chosen atoms."""
110 def __init__(self, indices=None, mask=None):
111 """Constrain chosen atoms.
113 Parameters
114 ----------
115 indices : list of int
116 Indices for those atoms that should be constrained.
117 mask : list of bool
118 One boolean per atom indicating if the atom should be
119 constrained or not.
121 Examples
122 --------
123 Fix all Copper atoms:
125 >>> mask = [s == 'Cu' for s in atoms.get_chemical_symbols()]
126 >>> c = FixAtoms(mask=mask)
127 >>> atoms.set_constraint(c)
129 Fix all atoms with z-coordinate less than 1.0 Angstrom:
131 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)
132 >>> atoms.set_constraint(c)
133 """
135 if indices is None and mask is None:
136 raise ValueError('Use "indices" or "mask".')
137 if indices is not None and mask is not None:
138 raise ValueError('Use only one of "indices" and "mask".')
140 if mask is not None:
141 indices = np.arange(len(mask))[np.asarray(mask, bool)]
142 else:
143 # Check for duplicates:
144 srt = np.sort(indices)
145 if (np.diff(srt) == 0).any():
146 raise ValueError(
147 'FixAtoms: The indices array contained duplicates. '
148 'Perhaps you wanted to specify a mask instead, but '
149 'forgot the mask= keyword.')
150 self.index = np.asarray(indices, int)
152 if self.index.ndim != 1:
153 raise ValueError('Wrong argument to FixAtoms class!')
155 def get_removed_dof(self, atoms):
156 return 3 * len(self.index)
158 def adjust_positions(self, atoms, new):
159 new[self.index] = atoms.positions[self.index]
161 def adjust_forces(self, atoms, forces):
162 forces[self.index] = 0.0
164 def index_shuffle(self, atoms, ind):
165 # See docstring of superclass
166 index = []
167 for new, old in slice2enlist(ind, len(atoms)):
168 if old in self.index:
169 index.append(new)
170 if len(index) == 0:
171 raise IndexError('All indices in FixAtoms not part of slice')
172 self.index = np.asarray(index, int)
174 def get_indices(self):
175 return self.index
177 def __repr__(self):
178 return 'FixAtoms(indices=%s)' % ints2string(self.index)
180 def todict(self):
181 return {'name': 'FixAtoms',
182 'kwargs': {'indices': self.index.tolist()}}
184 def repeat(self, m, n):
185 i0 = 0
186 natoms = 0
187 if isinstance(m, int):
188 m = (m, m, m)
189 index_new = []
190 for m2 in range(m[2]):
191 for m1 in range(m[1]):
192 for m0 in range(m[0]):
193 i1 = i0 + n
194 index_new += [i + natoms for i in self.index]
195 i0 = i1
196 natoms += n
197 self.index = np.asarray(index_new, int)
198 return self
200 def delete_atoms(self, indices, natoms):
201 """Removes atom number ind from the index array, if present.
203 Required for removing atoms with existing FixAtoms constraints.
204 """
206 i = np.zeros(natoms, int) - 1
207 new = np.delete(np.arange(natoms), indices)
208 i[new] = np.arange(len(new))
209 index = i[self.index]
210 self.index = index[index >= 0]
211 if len(self.index) == 0:
212 return None
213 return self
216class FixCom(FixConstraint):
217 """Constraint class for fixing the center of mass.
219 References
221 https://pubs.acs.org/doi/abs/10.1021/jp9722824
223 """
225 def get_removed_dof(self, atoms):
226 return 3
228 def adjust_positions(self, atoms, new):
229 masses = atoms.get_masses()
230 old_cm = atoms.get_center_of_mass()
231 new_cm = np.dot(masses, new) / masses.sum()
232 d = old_cm - new_cm
233 new += d
235 def adjust_forces(self, atoms, forces):
236 m = atoms.get_masses()
237 mm = np.tile(m, (3, 1)).T
238 lb = np.sum(mm * forces, axis=0) / sum(m**2)
239 forces -= mm * lb
241 def todict(self):
242 return {'name': 'FixCom',
243 'kwargs': {}}
246def ints2string(x, threshold=None):
247 """Convert ndarray of ints to string."""
248 if threshold is None or len(x) <= threshold:
249 return str(x.tolist())
250 return str(x[:threshold].tolist())[:-1] + ', ...]'
253class FixBondLengths(FixConstraint):
254 maxiter = 500
256 def __init__(self, pairs, tolerance=1e-13,
257 bondlengths=None, iterations=None):
258 """iterations:
259 Ignored"""
260 self.pairs = np.asarray(pairs)
261 self.tolerance = tolerance
262 self.bondlengths = bondlengths
264 def get_removed_dof(self, atoms):
265 return len(self.pairs)
267 def adjust_positions(self, atoms, new):
268 old = atoms.positions
269 masses = atoms.get_masses()
271 if self.bondlengths is None:
272 self.bondlengths = self.initialize_bond_lengths(atoms)
274 for i in range(self.maxiter):
275 converged = True
276 for j, ab in enumerate(self.pairs):
277 a = ab[0]
278 b = ab[1]
279 cd = self.bondlengths[j]
280 r0 = old[a] - old[b]
281 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
282 d1 = new[a] - new[b] - r0 + d0
283 m = 1 / (1 / masses[a] + 1 / masses[b])
284 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)
285 if abs(x) > self.tolerance:
286 new[a] += x * m / masses[a] * d0
287 new[b] -= x * m / masses[b] * d0
288 converged = False
289 if converged:
290 break
291 else:
292 raise RuntimeError('Did not converge')
294 def adjust_momenta(self, atoms, p):
295 old = atoms.positions
296 masses = atoms.get_masses()
298 if self.bondlengths is None:
299 self.bondlengths = self.initialize_bond_lengths(atoms)
301 for i in range(self.maxiter):
302 converged = True
303 for j, ab in enumerate(self.pairs):
304 a = ab[0]
305 b = ab[1]
306 cd = self.bondlengths[j]
307 d = old[a] - old[b]
308 d, _ = find_mic(d, atoms.cell, atoms.pbc)
309 dv = p[a] / masses[a] - p[b] / masses[b]
310 m = 1 / (1 / masses[a] + 1 / masses[b])
311 x = -np.dot(dv, d) / cd**2
312 if abs(x) > self.tolerance:
313 p[a] += x * m * d
314 p[b] -= x * m * d
315 converged = False
316 if converged:
317 break
318 else:
319 raise RuntimeError('Did not converge')
321 def adjust_forces(self, atoms, forces):
322 self.constraint_forces = -forces
323 self.adjust_momenta(atoms, forces)
324 self.constraint_forces += forces
326 def initialize_bond_lengths(self, atoms):
327 bondlengths = np.zeros(len(self.pairs))
329 for i, ab in enumerate(self.pairs):
330 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)
332 return bondlengths
334 def get_indices(self):
335 return np.unique(self.pairs.ravel())
337 def todict(self):
338 return {'name': 'FixBondLengths',
339 'kwargs': {'pairs': self.pairs.tolist(),
340 'tolerance': self.tolerance}}
342 def index_shuffle(self, atoms, ind):
343 """Shuffle the indices of the two atoms in this constraint"""
344 map = np.zeros(len(atoms), int)
345 map[ind] = 1
346 n = map.sum()
347 map[:] = -1
348 map[ind] = range(n)
349 pairs = map[self.pairs]
350 self.pairs = pairs[(pairs != -1).all(1)]
351 if len(self.pairs) == 0:
352 raise IndexError('Constraint not part of slice')
355def FixBondLength(a1, a2):
356 """Fix distance between atoms with indices a1 and a2."""
357 return FixBondLengths([(a1, a2)])
360class FixLinearTriatomic(FixConstraint):
361 """Holonomic constraints for rigid linear triatomic molecules."""
363 def __init__(self, triples):
364 """Apply RATTLE-type bond constraints between outer atoms n and m
365 and linear vectorial constraints to the position of central
366 atoms o to fix the geometry of linear triatomic molecules of the
367 type:
369 n--o--m
371 Parameters:
373 triples: list
374 Indices of the atoms forming the linear molecules to constrain
375 as triples. Sequence should be (n, o, m) or (m, o, n).
377 When using these constraints in molecular dynamics or structure
378 optimizations, atomic forces need to be redistributed within a
379 triple. The function redistribute_forces_optimization implements
380 the redistribution of forces for structure optimization, while
381 the function redistribute_forces_md implements the redistribution
382 for molecular dynamics.
384 References:
386 Ciccotti et al. Molecular Physics 47 (1982)
387 https://doi.org/10.1080/00268978200100942
388 """
389 self.triples = np.asarray(triples)
390 if self.triples.shape[1] != 3:
391 raise ValueError('"triples" has wrong size')
392 self.bondlengths = None
394 def get_removed_dof(self, atoms):
395 return 4 * len(self.triples)
397 @property
398 def n_ind(self):
399 return self.triples[:, 0]
401 @property
402 def m_ind(self):
403 return self.triples[:, 2]
405 @property
406 def o_ind(self):
407 return self.triples[:, 1]
409 def initialize(self, atoms):
410 masses = atoms.get_masses()
411 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)
413 self.bondlengths = self.initialize_bond_lengths(atoms)
414 self.bondlengths_nm = self.bondlengths.sum(axis=1)
416 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]
417 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +
418 C1[:, 1] ** 2 * self.mass_n * self.mass_o +
419 self.mass_n * self.mass_m)
420 C2 = C1 / C2[:, None]
421 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]
422 C3 = C2 * self.mass_o[:, None] * C3[:, None]
423 C3[:, 1] *= -1
424 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T
425 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)
426 C4 = C1 / C4[:, None]
428 self.C1 = C1
429 self.C2 = C2
430 self.C3 = C3
431 self.C4 = C4
433 def adjust_positions(self, atoms, new):
434 old = atoms.positions
435 new_n, new_m, new_o = self.get_slices(new)
437 if self.bondlengths is None:
438 self.initialize(atoms)
440 r0 = old[self.n_ind] - old[self.m_ind]
441 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)
442 d1 = new_n - new_m - r0 + d0
443 a = np.einsum('ij,ij->i', d0, d0)
444 b = np.einsum('ij,ij->i', d1, d0)
445 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2
446 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))
447 g = g[:, None] * self.C3
448 new_n -= g[:, 0, None] * d0
449 new_m += g[:, 1, None] * d0
450 if np.allclose(d0, r0):
451 new_o = (self.C1[:, 0, None] * new_n
452 + self.C1[:, 1, None] * new_m)
453 else:
454 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)
455 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)
456 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2
457 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)
459 self.set_slices(new_n, new_m, new_o, new)
461 def adjust_momenta(self, atoms, p):
462 old = atoms.positions
463 p_n, p_m, p_o = self.get_slices(p)
465 if self.bondlengths is None:
466 self.initialize(atoms)
468 mass_nn = self.mass_n[:, None]
469 mass_mm = self.mass_m[:, None]
470 mass_oo = self.mass_o[:, None]
472 d = old[self.n_ind] - old[self.m_ind]
473 d, _ = find_mic(d, atoms.cell, atoms.pbc)
474 dv = p_n / mass_nn - p_m / mass_mm
475 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2
476 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]
477 p_n -= k[:, 0, None] * mass_nn * d
478 p_m += k[:, 1, None] * mass_mm * d
479 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +
480 self.C1[:, 1, None] * p_m / mass_mm))
482 self.set_slices(p_n, p_m, p_o, p)
484 def adjust_forces(self, atoms, forces):
486 if self.bondlengths is None:
487 self.initialize(atoms)
489 A = self.C4 * np.diff(self.C1)
490 A[:, 0] *= -1
491 A -= 1
492 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]
493 A /= (A.sum(axis=1))[:, None]
495 self.constraint_forces = -forces
496 old = atoms.positions
498 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)
500 d = old[self.n_ind] - old[self.m_ind]
501 d, _ = find_mic(d, atoms.cell, atoms.pbc)
502 df = fr_n - fr_m
503 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2
504 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]
505 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]
506 forces[self.o_ind] = fr_o + k[:, None] * d * B
508 self.constraint_forces += forces
510 def redistribute_forces_optimization(self, forces):
511 """Redistribute forces within a triple when performing structure
512 optimizations.
514 The redistributed forces needs to be further adjusted using the
515 appropriate Lagrange multipliers as implemented in adjust_forces."""
516 forces_n, forces_m, forces_o = self.get_slices(forces)
517 C1_1 = self.C1[:, 0, None]
518 C1_2 = self.C1[:, 1, None]
519 C4_1 = self.C4[:, 0, None]
520 C4_2 = self.C4[:, 1, None]
522 fr_n = ((1 - C4_1 * C1_1) * forces_n -
523 C4_1 * (C1_2 * forces_m - forces_o))
524 fr_m = ((1 - C4_2 * C1_2) * forces_m -
525 C4_2 * (C1_1 * forces_n - forces_o))
526 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +
527 C4_1 * forces_n + C4_2 * forces_m)
529 return fr_n, fr_m, fr_o
531 def redistribute_forces_md(self, atoms, forces, rand=False):
532 """Redistribute forces within a triple when performing molecular
533 dynamics.
535 When rand=True, use the equations for random force terms, as
536 used e.g. by Langevin dynamics, otherwise apply the standard
537 equations for deterministic forces (see Ciccotti et al. Molecular
538 Physics 47 (1982))."""
539 if self.bondlengths is None:
540 self.initialize(atoms)
541 forces_n, forces_m, forces_o = self.get_slices(forces)
542 C1_1 = self.C1[:, 0, None]
543 C1_2 = self.C1[:, 1, None]
544 C2_1 = self.C2[:, 0, None]
545 C2_2 = self.C2[:, 1, None]
546 mass_nn = self.mass_n[:, None]
547 mass_mm = self.mass_m[:, None]
548 mass_oo = self.mass_o[:, None]
549 if rand:
550 mr1 = (mass_mm / mass_nn) ** 0.5
551 mr2 = (mass_oo / mass_nn) ** 0.5
552 mr3 = (mass_nn / mass_mm) ** 0.5
553 mr4 = (mass_oo / mass_mm) ** 0.5
554 else:
555 mr1 = 1.0
556 mr2 = 1.0
557 mr3 = 1.0
558 mr4 = 1.0
560 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -
561 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -
562 mr2 * mass_mm * mass_nn * forces_o))
564 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -
565 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -
566 mr4 * mass_mm * mass_nn * forces_o))
568 self.set_slices(fr_n, fr_m, 0.0, forces)
570 def get_slices(self, a):
571 a_n = a[self.n_ind]
572 a_m = a[self.m_ind]
573 a_o = a[self.o_ind]
575 return a_n, a_m, a_o
577 def set_slices(self, a_n, a_m, a_o, a):
578 a[self.n_ind] = a_n
579 a[self.m_ind] = a_m
580 a[self.o_ind] = a_o
582 def initialize_bond_lengths(self, atoms):
583 bondlengths = np.zeros((len(self.triples), 2))
585 for i in range(len(self.triples)):
586 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],
587 self.o_ind[i], mic=True)
588 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],
589 self.m_ind[i], mic=True)
591 return bondlengths
593 def get_indices(self):
594 return np.unique(self.triples.ravel())
596 def todict(self):
597 return {'name': 'FixLinearTriatomic',
598 'kwargs': {'triples': self.triples.tolist()}}
600 def index_shuffle(self, atoms, ind):
601 """Shuffle the indices of the three atoms in this constraint"""
602 map = np.zeros(len(atoms), int)
603 map[ind] = 1
604 n = map.sum()
605 map[:] = -1
606 map[ind] = range(n)
607 triples = map[self.triples]
608 self.triples = triples[(triples != -1).all(1)]
609 if len(self.triples) == 0:
610 raise IndexError('Constraint not part of slice')
613class FixedMode(FixConstraint):
614 """Constrain atoms to move along directions orthogonal to
615 a given mode only."""
617 def __init__(self, mode):
618 self.mode = (np.asarray(mode) / np.sqrt((mode**2).sum())).reshape(-1)
620 def get_removed_dof(self, atoms):
621 return len(atoms)
623 def adjust_positions(self, atoms, newpositions):
624 newpositions = newpositions.ravel()
625 oldpositions = atoms.positions.ravel()
626 step = newpositions - oldpositions
627 newpositions -= self.mode * np.dot(step, self.mode)
629 def adjust_forces(self, atoms, forces):
630 forces = forces.ravel()
631 forces -= self.mode * np.dot(forces, self.mode)
633 def index_shuffle(self, atoms, ind):
634 eps = 1e-12
635 mode = self.mode.reshape(-1, 3)
636 excluded = np.ones(len(mode), dtype=bool)
637 excluded[ind] = False
638 if (abs(mode[excluded]) > eps).any():
639 raise IndexError('All nonzero parts of mode not in slice')
640 self.mode = mode[ind].ravel()
642 def get_indices(self):
643 # This function will never properly work because it works on all
644 # atoms and it has no idea how to tell how many atoms it is
645 # attached to. If it is being used, surely the user knows
646 # everything is being constrained.
647 return []
649 def todict(self):
650 return {'name': 'FixedMode',
651 'kwargs': {'mode': self.mode.tolist()}}
653 def __repr__(self):
654 return 'FixedMode(%s)' % self.mode.tolist()
657class FixedPlane(FixConstraintSingle):
658 """Constrain an atom index *a* to move in a given plane only.
660 The plane is defined by its normal vector *direction*."""
662 def __init__(self, a, direction):
663 self.a = a
664 self.dir = np.asarray(direction) / sqrt(np.dot(direction, direction))
666 def get_removed_dof(self, atoms):
667 return 1
669 def adjust_positions(self, atoms, newpositions):
670 step = newpositions[self.a] - atoms.positions[self.a]
671 newpositions[self.a] -= self.dir * np.dot(step, self.dir)
673 def adjust_forces(self, atoms, forces):
674 forces[self.a] -= self.dir * np.dot(forces[self.a], self.dir)
676 def todict(self):
677 return {'name': 'FixedPlane',
678 'kwargs': {'a': self.a, 'direction': self.dir.tolist()}}
680 def __repr__(self):
681 return 'FixedPlane(%d, %s)' % (self.a, self.dir.tolist())
684class FixedLine(FixConstraintSingle):
685 """Constrain an atom index *a* to move on a given line only.
687 The line is defined by its vector *direction*."""
689 def __init__(self, a, direction):
690 self.a = a
691 self.dir = np.asarray(direction) / sqrt(np.dot(direction, direction))
693 def get_removed_dof(self, atoms):
694 return 2
696 def adjust_positions(self, atoms, newpositions):
697 step = newpositions[self.a] - atoms.positions[self.a]
698 x = np.dot(step, self.dir)
699 newpositions[self.a] = atoms.positions[self.a] + x * self.dir
701 def adjust_forces(self, atoms, forces):
702 forces[self.a] = self.dir * np.dot(forces[self.a], self.dir)
704 def __repr__(self):
705 return 'FixedLine(%d, %s)' % (self.a, self.dir.tolist())
707 def todict(self):
708 return {'name': 'FixedLine',
709 'kwargs': {'a': self.a, 'direction': self.dir.tolist()}}
712class FixCartesian(FixConstraintSingle):
713 'Fix an atom index *a* in the directions of the cartesian coordinates.'
715 def __init__(self, a, mask=(1, 1, 1)):
716 self.a = a
717 self.mask = ~np.asarray(mask, bool)
719 def get_removed_dof(self, atoms):
720 return 3 - self.mask.sum()
722 def adjust_positions(self, atoms, new):
723 step = new[self.a] - atoms.positions[self.a]
724 step *= self.mask
725 new[self.a] = atoms.positions[self.a] + step
727 def adjust_forces(self, atoms, forces):
728 forces[self.a] *= self.mask
730 def __repr__(self):
731 return 'FixCartesian(a={0}, mask={1})'.format(self.a,
732 list(~self.mask))
734 def todict(self):
735 return {'name': 'FixCartesian',
736 'kwargs': {'a': self.a, 'mask': ~self.mask.tolist()}}
739class FixScaled(FixConstraintSingle):
740 'Fix an atom index *a* in the directions of the unit vectors.'
742 def __init__(self, cell, a, mask=(1, 1, 1)):
743 self.cell = np.asarray(cell)
744 self.a = a
745 self.mask = np.array(mask, bool)
747 def get_removed_dof(self, atoms):
748 return self.mask.sum()
750 def adjust_positions(self, atoms, new):
751 scaled_old = atoms.cell.scaled_positions(atoms.positions)
752 scaled_new = atoms.cell.scaled_positions(new)
753 for n in range(3):
754 if self.mask[n]:
755 scaled_new[self.a, n] = scaled_old[self.a, n]
756 new[self.a] = atoms.cell.cartesian_positions(scaled_new)[self.a]
758 def adjust_forces(self, atoms, forces):
759 # Forces are covarient to the coordinate transformation,
760 # use the inverse transformations
761 scaled_forces = atoms.cell.cartesian_positions(forces)
762 scaled_forces[self.a] *= -(self.mask - 1)
763 forces[self.a] = atoms.cell.scaled_positions(scaled_forces)[self.a]
765 def todict(self):
766 return {'name': 'FixScaled',
767 'kwargs': {'a': self.a,
768 'cell': self.cell.tolist(),
769 'mask': self.mask.tolist()}}
771 def __repr__(self):
772 return 'FixScaled(%s, %d, %s)' % (repr(self.cell),
773 self.a,
774 repr(self.mask))
777# TODO: Better interface might be to use dictionaries in place of very
778# nested lists/tuples
779class FixInternals(FixConstraint):
780 """Constraint object for fixing multiple internal coordinates.
782 Allows fixing bonds, angles, and dihedrals.
783 Please provide angular units in degrees using angles_deg and
784 dihedrals_deg.
785 """
786 def __init__(self, bonds=None, angles=None, dihedrals=None,
787 angles_deg=None, dihedrals_deg=None,
788 bondcombos=None, anglecombos=None, dihedralcombos=None,
789 mic=False, epsilon=1.e-7):
791 # deprecate public API using radians; degrees is preferred
792 warn_msg = 'Please specify {} in degrees using the {} argument.'
793 if angles:
794 warn(FutureWarning(warn_msg.format('angles', 'angle_deg')))
795 angles = np.asarray(angles)
796 angles[:, 0] = angles[:, 0] / np.pi * 180
797 angles = angles.tolist()
798 else:
799 angles = angles_deg
800 if dihedrals:
801 warn(FutureWarning(warn_msg.format('dihedrals', 'dihedrals_deg')))
802 dihedrals = np.asarray(dihedrals)
803 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180
804 dihedrals = dihedrals.tolist()
805 else:
806 dihedrals = dihedrals_deg
808 self.bonds = bonds or []
809 self.angles = angles or []
810 self.dihedrals = dihedrals or []
811 self.bondcombos = bondcombos or []
812 self.anglecombos = anglecombos or []
813 self.dihedralcombos = dihedralcombos or []
814 self.mic = mic
815 self.epsilon = epsilon
817 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)
818 + len(self.bondcombos) + len(self.anglecombos)
819 + len(self.dihedralcombos))
821 # Initialize these at run-time:
822 self.constraints = []
823 self.initialized = False
825 def get_removed_dof(self, atoms):
826 return self.n
828 def initialize(self, atoms):
829 if self.initialized:
830 return
831 masses = np.repeat(atoms.get_masses(), 3)
832 cell = None
833 pbc = None
834 if self.mic:
835 cell = atoms.cell
836 pbc = atoms.pbc
837 self.constraints = []
838 for data, make_constr in [(self.bonds, self.FixBondLengthAlt),
839 (self.angles, self.FixAngle),
840 (self.dihedrals, self.FixDihedral),
841 (self.bondcombos, self.FixBondCombo),
842 (self.anglecombos, self.FixAngleCombo),
843 (self.dihedralcombos, self.FixDihedralCombo)]:
844 for datum in data:
845 constr = make_constr(datum[0], datum[1], masses, cell, pbc)
846 self.constraints.append(constr)
847 self.initialized = True
849 def shuffle_definitions(self, shuffle_dic, internal_type):
850 dfns = [] # definitions
851 for dfn in internal_type: # e.g. for bond in self.bonds
852 append = True
853 new_dfn = [dfn[0], list(dfn[1])]
854 for old in dfn[1]:
855 if old in shuffle_dic:
856 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]
857 else:
858 append = False
859 break
860 if append:
861 dfns.append(new_dfn)
862 return dfns
864 def shuffle_combos(self, shuffle_dic, internal_type):
865 dfns = [] # definitions
866 for dfn in internal_type: # e.g. for bondcombo in self.bondcombos
867 append = True
868 all_indices = [idx[0:-1] for idx in dfn[1]]
869 new_dfn = [dfn[0], list(dfn[1])]
870 for i, indices in enumerate(all_indices):
871 for old in indices:
872 if old in shuffle_dic:
873 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]
874 else:
875 append = False
876 break
877 if not append:
878 break
879 if append:
880 dfns.append(new_dfn)
881 return dfns
883 def index_shuffle(self, atoms, ind):
884 # See docstring of superclass
885 self.initialize(atoms)
886 shuffle_dic = dict(slice2enlist(ind, len(atoms)))
887 shuffle_dic = {old: new for new, old in shuffle_dic.items()}
888 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)
889 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)
890 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)
891 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)
892 self.anglecombos = self.shuffle_combos(shuffle_dic, self.anglecombos)
893 self.dihedralcombos = self.shuffle_combos(shuffle_dic,
894 self.dihedralcombos)
895 self.initialized = False
896 self.initialize(atoms)
897 if len(self.constraints) == 0:
898 raise IndexError('Constraint not part of slice')
900 def get_indices(self):
901 cons = []
902 for dfn in self.bonds + self.dihedrals + self.angles:
903 cons.extend(dfn[1])
904 for dfn in self.bondcombos + self.anglecombos + self.dihedralcombos:
905 for partial_dfn in dfn[1]:
906 cons.extend(partial_dfn[0:-1]) # last index is the coefficient
907 return list(set(cons))
909 def todict(self):
910 return {'name': 'FixInternals',
911 'kwargs': {'bonds': self.bonds,
912 'angles': self.angles,
913 'dihedrals': self.dihedrals,
914 'bondcombos': self.bondcombos,
915 'anglecombos': self.anglecombos,
916 'dihedralcombos': self.dihedralcombos,
917 'mic': self.mic,
918 'epsilon': self.epsilon}}
920 def adjust_positions(self, atoms, new):
921 self.initialize(atoms)
922 for constraint in self.constraints:
923 constraint.prepare_jacobian(atoms.positions)
924 for j in range(50):
925 maxerr = 0.0
926 for constraint in self.constraints:
927 constraint.adjust_positions(atoms.positions, new)
928 maxerr = max(abs(constraint.sigma), maxerr)
929 if maxerr < self.epsilon:
930 return
931 raise ValueError('Shake did not converge.')
933 def adjust_forces(self, atoms, forces):
934 """Project out translations and rotations and all other constraints"""
935 self.initialize(atoms)
936 positions = atoms.positions
937 N = len(forces)
938 list2_constraints = list(np.zeros((6, N, 3)))
939 tx, ty, tz, rx, ry, rz = list2_constraints
941 list_constraints = [r.ravel() for r in list2_constraints]
943 tx[:, 0] = 1.0
944 ty[:, 1] = 1.0
945 tz[:, 2] = 1.0
946 ff = forces.ravel()
948 # Calculate the center of mass
949 center = positions.sum(axis=0) / N
951 rx[:, 1] = -(positions[:, 2] - center[2])
952 rx[:, 2] = positions[:, 1] - center[1]
953 ry[:, 0] = positions[:, 2] - center[2]
954 ry[:, 2] = -(positions[:, 0] - center[0])
955 rz[:, 0] = -(positions[:, 1] - center[1])
956 rz[:, 1] = positions[:, 0] - center[0]
958 # Normalizing transl., rotat. constraints
959 for r in list2_constraints:
960 r /= np.linalg.norm(r.ravel())
962 # Add all angle, etc. constraint vectors
963 for constraint in self.constraints:
964 constraint.prepare_jacobian(positions)
965 constraint.adjust_forces(positions, forces)
966 list_constraints.insert(0, constraint.jacobian)
967 # QR DECOMPOSITION - GRAM SCHMIDT
969 list_constraints = [r.ravel() for r in list_constraints]
970 aa = np.column_stack(list_constraints)
971 (aa, bb) = np.linalg.qr(aa)
972 # Projection
973 hh = []
974 for i, constraint in enumerate(self.constraints):
975 hh.append(aa[:, i] * np.row_stack(aa[:, i]))
977 txx = aa[:, self.n] * np.row_stack(aa[:, self.n])
978 tyy = aa[:, self.n + 1] * np.row_stack(aa[:, self.n + 1])
979 tzz = aa[:, self.n + 2] * np.row_stack(aa[:, self.n + 2])
980 rxx = aa[:, self.n + 3] * np.row_stack(aa[:, self.n + 3])
981 ryy = aa[:, self.n + 4] * np.row_stack(aa[:, self.n + 4])
982 rzz = aa[:, self.n + 5] * np.row_stack(aa[:, self.n + 5])
983 T = txx + tyy + tzz + rxx + ryy + rzz
984 for vec in hh:
985 T += vec
986 ff = np.dot(T, np.row_stack(ff))
987 forces[:, :] -= np.dot(T, np.row_stack(ff)).reshape(-1, 3)
989 def __repr__(self):
990 constraints = repr(self.constraints)
991 return 'FixInternals(_copy_init=%s, epsilon=%s)' % (constraints,
992 repr(self.epsilon))
994 def __str__(self):
995 return '\n'.join([repr(c) for c in self.constraints])
997 # Classes for internal use in FixInternals
998 class FixInternalsBase:
999 """Base class for subclasses of FixInternals."""
1000 def __init__(self, targetvalue, indices, masses, cell, pbc):
1001 self.targetvalue = targetvalue # constant target value
1002 self.indices = [defin[0:-1] for defin in indices] # indices, defs
1003 self.coefs = np.asarray([defin[-1] for defin in indices]) # coefs
1004 self.masses = masses
1005 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix
1006 self.sigma = 1. # difference between current and target value
1007 self.projected_force = None # helps optimizers scan along constr.
1008 self.cell = cell
1009 self.pbc = pbc
1011 def finalize_jacobian(self, pos, n_internals, n, derivs):
1012 """Populate jacobian with derivatives for `n_internals` defined
1013 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""
1014 jacobian = np.zeros((n_internals, *pos.shape))
1015 for i, idx in enumerate(self.indices):
1016 for j in range(n):
1017 jacobian[i, idx[j]] = derivs[i, j]
1018 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))
1019 self.jacobian = self.coefs @ jacobian
1021 def finalize_positions(self, newpos):
1022 jacobian = self.jacobian / self.masses
1023 lamda = -self.sigma / np.dot(jacobian, self.jacobian)
1024 dnewpos = lamda * jacobian
1025 newpos += dnewpos.reshape(newpos.shape)
1027 def adjust_forces(self, positions, forces):
1028 self.projected_force = np.dot(self.jacobian, forces.ravel())
1029 self.jacobian /= np.linalg.norm(self.jacobian)
1031 class FixBondCombo(FixInternalsBase):
1032 """Constraint subobject for fixing linear combination of bond lengths
1033 within FixInternals.
1035 sum_i( coef_i * bond_length_i ) = constant
1036 """
1037 def prepare_jacobian(self, pos):
1038 bondvectors = [pos[k] - pos[h] for h, k in self.indices]
1039 derivs = get_distances_derivatives(bondvectors, cell=self.cell,
1040 pbc=self.pbc)
1041 self.finalize_jacobian(pos, len(bondvectors), 2, derivs)
1043 def adjust_positions(self, oldpos, newpos):
1044 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]
1045 (_, ), (dists, ) = conditional_find_mic([bondvectors],
1046 cell=self.cell,
1047 pbc=self.pbc)
1048 value = np.dot(self.coefs, dists)
1049 self.sigma = value - self.targetvalue
1050 self.finalize_positions(newpos)
1052 def __repr__(self):
1053 return 'FixBondCombo({}, {}, {})'.format(repr(self.targetvalue),
1054 self.indices, self.coefs)
1056 class FixBondLengthAlt(FixBondCombo):
1057 """Constraint subobject for fixing bond length within FixInternals.
1058 Fix distance between atoms with indices a1, a2."""
1059 def __init__(self, targetvalue, indices, masses, cell, pbc):
1060 indices = [list(indices) + [1.]] # bond definition with coef 1.
1061 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1063 def __repr__(self):
1064 return 'FixBondLengthAlt({}, {})'.format(self.targetvalue,
1065 *self.indices)
1067 class FixAngleCombo(FixInternalsBase):
1068 """Constraint subobject for fixing linear combination of angles
1069 within FixInternals.
1071 sum_i( coef_i * angle_i ) = constant
1072 """
1073 def gather_vectors(self, pos):
1074 v0 = [pos[h] - pos[k] for h, k, l in self.indices]
1075 v1 = [pos[l] - pos[k] for h, k, l in self.indices]
1076 return v0, v1
1078 def prepare_jacobian(self, pos):
1079 v0, v1 = self.gather_vectors(pos)
1080 derivs = get_angles_derivatives(v0, v1, cell=self.cell,
1081 pbc=self.pbc)
1082 self.finalize_jacobian(pos, len(v0), 3, derivs)
1084 def adjust_positions(self, oldpos, newpos):
1085 v0, v1 = self.gather_vectors(newpos)
1086 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)
1087 value = np.dot(self.coefs, value)
1088 self.sigma = value - self.targetvalue
1089 self.finalize_positions(newpos)
1091 def __repr__(self):
1092 return 'FixAngleCombo({}, {}, {})'.format(self.targetvalue,
1093 self.indices, self.coefs)
1095 class FixAngle(FixAngleCombo):
1096 """Constraint object for fixing an angle within
1097 FixInternals using the SHAKE algorithm.
1099 SHAKE convergence is potentially problematic for angles very close to
1100 0 or 180 degrees as there is a singularity in the Cartesian derivative.
1101 """
1102 def __init__(self, targetvalue, indices, masses, cell, pbc):
1103 """Fix atom movement to construct a constant angle."""
1104 indices = [list(indices) + [1.]] # angle definition with coef 1.
1105 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1107 def __repr__(self):
1108 return 'FixAngle({}, {})'.format(self.targetvalue, *self.indices)
1110 class FixDihedralCombo(FixInternalsBase):
1111 """Constraint subobject for fixing linear combination of dihedrals
1112 within FixInternals.
1114 sum_i( coef_i * dihedral_i ) = constant
1115 """
1116 def gather_vectors(self, pos):
1117 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]
1118 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]
1119 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]
1120 return v0, v1, v2
1122 def prepare_jacobian(self, pos):
1123 v0, v1, v2 = self.gather_vectors(pos)
1124 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,
1125 pbc=self.pbc)
1126 self.finalize_jacobian(pos, len(v0), 4, derivs)
1128 def adjust_positions(self, oldpos, newpos):
1129 v0, v1, v2 = self.gather_vectors(newpos)
1130 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1131 value = np.dot(self.coefs, value)
1132 self.sigma = value - self.targetvalue
1133 self.finalize_positions(newpos)
1135 def __repr__(self):
1136 return 'FixDihedralCombo({}, {}, {})'.format(self.targetvalue,
1137 self.indices,
1138 self.coefs)
1140 class FixDihedral(FixDihedralCombo):
1141 """Constraint object for fixing a dihedral angle using
1142 the SHAKE algorithm. This one allows also other constraints.
1144 SHAKE convergence is potentially problematic for near-undefined
1145 dihedral angles (i.e. when one of the two angles a012 or a123
1146 approaches 0 or 180 degrees).
1147 """
1148 def __init__(self, targetvalue, indices, masses, cell, pbc):
1149 indices = [list(indices) + [1.]] # dihedral def. with coef 1.
1150 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)
1152 def adjust_positions(self, oldpos, newpos):
1153 v0, v1, v2 = self.gather_vectors(newpos)
1154 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)
1155 # apply minimum dihedral difference 'convention': (diff <= 180)
1156 self.sigma = (value - self.targetvalue + 180) % 360 - 180
1157 self.finalize_positions(newpos)
1159 def __repr__(self):
1160 return 'FixDihedral({}, {})'.format(self.targetvalue, *self.indices)
1163class FixParametricRelations(FixConstraint):
1165 def __init__(
1166 self,
1167 indices,
1168 Jacobian,
1169 const_shift,
1170 params=None,
1171 eps=1e-12,
1172 use_cell=False,
1173 ):
1174 """Constrains the degrees of freedom to act in a reduced parameter space defined by the Jacobian
1176 These constraints are based off the work in: https://arxiv.org/abs/1908.01610
1178 The constraints linearly maps the full 3N degrees of freedom, where N is number of active
1179 lattice vectors/atoms onto a reduced subset of M free parameters, where M <= 3*N. The
1180 Jacobian matrix and constant shift vector map the full set of degrees of freedom onto the
1181 reduced parameter space.
1183 Currently the constraint is set up to handle either atomic positions or lattice vectors
1184 at one time, but not both. To do both simply add a two constraints for each set. This is
1185 done to keep the mathematics behind the operations separate.
1187 It would be possible to extend these constraints to allow non-linear transformations
1188 if functionality to update the Jacobian at each position update was included. This would
1189 require passing an update function evaluate it every time adjust_positions is callled.
1190 This is currently NOT supported, and there are no plans to implement it in the future.
1192 Args:
1193 indices (list of int): indices of the constrained atoms
1194 (if not None or empty then cell_indices must be None or Empty)
1195 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))): The Jacobian describing
1196 the parameter space transformation
1197 const_shift (np.ndarray(shape=(3*len(indices)))): A vector describing the constant term
1198 in the transformation not accounted for in the Jacobian
1199 params (list of str): parameters used in the parametric representation
1200 if None a list is generated based on the shape of the Jacobian
1201 eps (float): a small number to compare the similarity of numbers and set the precision used
1202 to generate the constraint expressions
1203 use_cell (bool): if True then act on the cell object
1204 """
1205 self.indices = np.array(indices)
1206 self.Jacobian = np.array(Jacobian)
1207 self.const_shift = np.array(const_shift)
1209 assert self.const_shift.shape[0] == 3*len(self.indices)
1210 assert self.Jacobian.shape[0] == 3*len(self.indices)
1212 self.eps = eps
1213 self.use_cell = use_cell
1215 if params is None:
1216 params = []
1217 if self.Jacobian.shape[1] > 0:
1218 int_fmt_str = "{:0" + str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"
1219 for param_ind in range(self.Jacobian.shape[1]):
1220 params.append("param_" + int_fmt_str.format(param_ind))
1221 else:
1222 assert len(params) == self.Jacobian.shape[-1]
1224 self.params = params
1226 self.Jacobian_inv = np.linalg.inv(self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T
1228 @classmethod
1229 def from_expressions(cls, indices, params, expressions, eps=1e-12, use_cell=False):
1230 """Converts the expressions into a Jacobian Matrix/const_shift vector and constructs a FixParametricRelations constraint
1232 The expressions must be a list like object of size 3*N and elements must be ordered as:
1233 [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],
1234 where i, j, and k are the first, second and third component of the atomic position/lattice
1235 vector. Currently only linear operations are allowed to be included in the expressions so
1236 only terms like:
1237 - const * param_0
1238 - sqrt[const] * param_1
1239 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M
1240 where const is any real number and param_0, param_1, ..., param_M are the parameters passed in
1241 params, are allowed.
1243 For example, the fractional atomic position constraints for wurtzite are:
1244 params = ["z1", "z2"]
1245 expressions = [
1246 "1.0/3.0", "2.0/3.0", "z1",
1247 "2.0/3.0", "1.0/3.0", "0.5 + z1",
1248 "1.0/3.0", "2.0/3.0", "z2",
1249 "2.0/3.0", "1.0/3.0", "0.5 + z2",
1250 ]
1252 For diamond are:
1253 params = []
1254 expressions = [
1255 "0.0", "0.0", "0.0",
1256 "0.25", "0.25", "0.25",
1257 ],
1259 and for stannite are
1260 params=["x4", "z4"]
1261 expressions = [
1262 "0.0", "0.0", "0.0",
1263 "0.0", "0.5", "0.5",
1264 "0.75", "0.25", "0.5",
1265 "0.25", "0.75", "0.5",
1266 "x4 + z4", "x4 + z4", "2*x4",
1267 "x4 - z4", "x4 - z4", "-2*x4",
1268 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",
1269 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",
1270 ]
1272 Args:
1273 indices (list of int): indices of the constrained atoms
1274 (if not None or empty then cell_indices must be None or Empty)
1275 params (list of str): parameters used in the parametric representation
1276 expressions (list of str): expressions used to convert from the parametric to the real space
1277 representation
1278 eps (float): a small number to compare the similarity of numbers and set the precision used
1279 to generate the constraint expressions
1280 use_cell (bool): if True then act on the cell object
1282 Returns:
1283 cls(
1284 indices,
1285 Jacobian generated from expressions,
1286 const_shift generated from expressions,
1287 params,
1288 eps-12,
1289 use_cell,
1290 )
1291 """
1292 Jacobian = np.zeros((3*len(indices), len(params)))
1293 const_shift = np.zeros(3*len(indices))
1295 for expr_ind, expression in enumerate(expressions):
1296 expression = expression.strip()
1298 # Convert subtraction to addition
1299 expression = expression.replace("-", "+(-1.0)*")
1300 if "+" == expression[0]:
1301 expression = expression[1:]
1302 elif "(+" == expression[:2]:
1303 expression = "(" + expression[2:]
1305 # Explicitly add leading zeros so when replacing param_1 with 0.0 param_11 does not become 0.01
1306 int_fmt_str = "{:0" + str(int(np.ceil(np.log10(len(params)+1)))) + "d}"
1308 param_dct = dict()
1309 param_map = dict()
1311 # Construct a standardized param template for A/B filling
1312 for param_ind, param in enumerate(params):
1313 param_str = "param_" + int_fmt_str.format(param_ind)
1314 param_map[param] = param_str
1315 param_dct[param_str] = 0.0
1317 # Replace the parameters according to the map
1318 # Sort by string length (long to short) to prevent cases like x11 becoming f"{param_map["x1"]}1"
1319 for param in sorted(params, key=lambda s: -1.0*len(s)):
1320 expression = expression.replace(param, param_map[param])
1322 # Partial linearity check
1323 for express_sec in expression.split("+"):
1324 in_sec = [param in express_sec for param in param_dct]
1325 n_params_in_sec = len(np.where(np.array(in_sec))[0])
1326 if n_params_in_sec > 1:
1327 raise ValueError("The FixParametricRelations expressions must be linear.")
1329 const_shift[expr_ind] = float(eval_expression(expression, param_dct))
1331 for param_ind in range(len(params)):
1332 param_str = "param_" + int_fmt_str.format(param_ind)
1333 if param_str not in expression:
1334 Jacobian[expr_ind, param_ind] = 0.0
1335 continue
1336 param_dct[param_str] = 1.0
1337 test_1 = float(eval_expression(expression, param_dct))
1338 test_1 -= const_shift[expr_ind]
1339 Jacobian[expr_ind, param_ind] = test_1
1341 param_dct[param_str] = 2.0
1342 test_2 = float(eval_expression(expression, param_dct))
1343 test_2 -= const_shift[expr_ind]
1344 if abs(test_2 / test_1 - 2.0) > eps:
1345 raise ValueError("The FixParametricRelations expressions must be linear.")
1346 param_dct[param_str] = 0.0
1348 args = [
1349 indices,
1350 Jacobian,
1351 const_shift,
1352 params,
1353 eps,
1354 use_cell,
1355 ]
1356 if cls is FixScaledParametricRelations:
1357 args = args[:-1]
1358 return cls(*args)
1360 @property
1361 def expressions(self):
1362 """Generate the expressions represented by the current self.Jacobian and self.const_shift objects"""
1363 expressions = []
1364 per = int(round(-1 * np.log10(self.eps)))
1365 fmt_str = "{:." + str(per + 1) + "g}"
1366 for index, shift_val in enumerate(self.const_shift):
1367 exp = ""
1368 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(shift_val) > self.eps:
1369 exp += fmt_str.format(shift_val)
1371 param_exp = ""
1372 for param_index, jacob_val in enumerate(self.Jacobian[index]):
1373 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)
1374 if abs_jacob_val < self.eps:
1375 continue
1377 param = self.params[param_index]
1378 if param_exp or exp:
1379 if jacob_val > -1.0*self.eps:
1380 param_exp += " + "
1381 else:
1382 param_exp += " - "
1383 elif (not exp) and (not param_exp) and (jacob_val < -1.0*self.eps):
1384 param_exp += "-"
1386 if np.abs(abs_jacob_val-1.0) <= self.eps:
1387 param_exp += "{:s}".format(param)
1388 else:
1389 param_exp += (fmt_str + "*{:s}").format(abs_jacob_val, param)
1391 exp += param_exp
1393 expressions.append(exp)
1394 return np.array(expressions).reshape((-1, 3))
1396 def todict(self):
1397 """Create a dictionary representation of the constraint"""
1398 return {
1399 "name": type(self).__name__,
1400 "kwargs": {
1401 "indices": self.indices,
1402 "params": self.params,
1403 "Jacobian": self.Jacobian,
1404 "const_shift": self.const_shift,
1405 "eps": self.eps,
1406 "use_cell": self.use_cell,
1407 }
1408 }
1410 def __repr__(self):
1411 """The str representation of the constraint"""
1412 if len(self.indices) > 1:
1413 indices_str = "[{:d}, ..., {:d}]".format(self.indices[0], self.indices[-1])
1414 else:
1415 indices_str = "[{:d}]".format(self.indices[0])
1417 if len(self.params) > 1:
1418 params_str = "[{:s}, ..., {:s}]".format(self.params[0], self.params[-1])
1419 elif len(self.params) == 1:
1420 params_str = "[{:s}]".format(self.params[0])
1421 else:
1422 params_str = "[]"
1424 return '{:s}({:s}, {:s}, ..., {:e})'.format(
1425 type(self).__name__,
1426 indices_str,
1427 params_str,
1428 self.eps
1429 )
1432class FixScaledParametricRelations(FixParametricRelations):
1434 def __init__(
1435 self,
1436 indices,
1437 Jacobian,
1438 const_shift,
1439 params=None,
1440 eps=1e-12,
1441 ):
1442 """The fractional coordinate version of FixParametricRelations
1444 All arguments are the same, but since this is for fractional coordinates use_cell is false
1445 """
1446 super(FixScaledParametricRelations, self).__init__(
1447 indices,
1448 Jacobian,
1449 const_shift,
1450 params,
1451 eps,
1452 False,
1453 )
1455 def adjust_contravariant(self, cell, vecs, B):
1456 """Adjust the values of a set of vectors that are contravariant with the unit transformation"""
1457 scaled = cell.scaled_positions(vecs).flatten()
1458 scaled = self.Jacobian_inv @ (scaled - B)
1459 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))
1461 return cell.cartesian_positions(scaled)
1463 def adjust_positions(self, atoms, positions):
1464 """Adjust positions of the atoms to match the constraints"""
1465 positions[self.indices] = self.adjust_contravariant(
1466 atoms.cell,
1467 positions[self.indices],
1468 self.const_shift,
1469 )
1470 positions[self.indices] = self.adjust_B(atoms.cell, positions[self.indices])
1472 def adjust_B(self, cell, positions):
1473 """Wraps the positions back to the unit cell and adjust B to keep track of this change"""
1474 fractional = cell.scaled_positions(positions)
1475 wrapped_fractional = (fractional % 1.0) % 1.0
1476 self.const_shift += np.round(wrapped_fractional - fractional).flatten()
1477 return cell.cartesian_positions(wrapped_fractional)
1479 def adjust_momenta(self, atoms, momenta):
1480 """Adjust momenta of the atoms to match the constraints"""
1481 momenta[self.indices] = self.adjust_contravariant(
1482 atoms.cell,
1483 momenta[self.indices],
1484 np.zeros(self.const_shift.shape),
1485 )
1487 def adjust_forces(self, atoms, forces):
1488 """Adjust forces of the atoms to match the constraints"""
1489 # Forces are coavarient to the coordinate transformation, use the inverse transformations
1490 cart2frac_jacob = np.zeros(2*(3*len(atoms),))
1491 for i_atom in range(len(atoms)):
1492 cart2frac_jacob[3*i_atom:3*(i_atom+1), 3*i_atom:3*(i_atom+1)] = atoms.cell.T
1494 jacobian = cart2frac_jacob @ self.Jacobian
1495 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T
1497 reduced_forces = jacobian.T @ forces.flatten()
1498 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)
1500 def todict(self):
1501 """Create a dictionary representation of the constraint"""
1502 dct = super(FixScaledParametricRelations, self).todict()
1503 del(dct["kwargs"]["use_cell"])
1504 return dct
1507class FixCartesianParametricRelations(FixParametricRelations):
1509 def __init__(
1510 self,
1511 indices,
1512 Jacobian,
1513 const_shift,
1514 params=None,
1515 eps=1e-12,
1516 use_cell=False,
1517 ):
1518 """The Cartesian coordinate version of FixParametricRelations"""
1519 super(FixCartesianParametricRelations, self).__init__(
1520 indices,
1521 Jacobian,
1522 const_shift,
1523 params,
1524 eps,
1525 use_cell,
1526 )
1528 def adjust_contravariant(self, vecs, B):
1529 """Adjust the values of a set of vectors that are contravariant with the unit transformation"""
1530 vecs = self.Jacobian_inv @ (vecs.flatten() - B)
1531 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))
1532 return vecs
1534 def adjust_positions(self, atoms, positions):
1535 """Adjust positions of the atoms to match the constraints"""
1536 if self.use_cell:
1537 return
1538 positions[self.indices] = self.adjust_contravariant(
1539 positions[self.indices],
1540 self.const_shift,
1541 )
1543 def adjust_momenta(self, atoms, momenta):
1544 """Adjust momenta of the atoms to match the constraints"""
1545 if self.use_cell:
1546 return
1547 momenta[self.indices] = self.adjust_contravariant(
1548 momenta[self.indices],
1549 np.zeros(self.const_shift.shape),
1550 )
1552 def adjust_forces(self, atoms, forces):
1553 """Adjust forces of the atoms to match the constraints"""
1554 if self.use_cell:
1555 return
1557 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()
1558 forces[self.indices] = (self.Jacobian_inv.T @ forces_reduced).reshape(-1, 3)
1560 def adjust_cell(self, atoms, cell):
1561 """Adjust the cell of the atoms to match the constraints"""
1562 if not self.use_cell:
1563 return
1564 cell[self.indices] = self.adjust_contravariant(
1565 cell[self.indices],
1566 np.zeros(self.const_shift.shape),
1567 )
1569 def adjust_stress(self, atoms, stress):
1570 """Adjust the stress of the atoms to match the constraints"""
1571 if not self.use_cell:
1572 return
1574 stress_3x3 = voigt_6_to_full_3x3_stress(stress)
1575 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()
1576 stress_3x3[self.indices] = (self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)
1578 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)
1581class Hookean(FixConstraint):
1582 """Applies a Hookean restorative force between a pair of atoms, an atom
1583 and a point, or an atom and a plane."""
1585 def __init__(self, a1, a2, k, rt=None):
1586 """Forces two atoms to stay close together by applying no force if
1587 they are below a threshold length, rt, and applying a Hookean
1588 restorative force when the distance between them exceeds rt. Can
1589 also be used to tether an atom to a fixed point in space or to a
1590 distance above a plane.
1592 a1 : int
1593 Index of atom 1
1594 a2 : one of three options
1595 1) index of atom 2
1596 2) a fixed point in cartesian space to which to tether a1
1597 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.
1598 k : float
1599 Hooke's law (spring) constant to apply when distance
1600 exceeds threshold_length. Units of eV A^-2.
1601 rt : float
1602 The threshold length below which there is no force. The
1603 length is 1) between two atoms, 2) between atom and point.
1604 This argument is not supplied in case 3. Units of A.
1606 If a plane is specified, the Hooke's law force is applied if the atom
1607 is on the normal side of the plane. For instance, the plane with
1608 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z
1609 intercept of +7 and a normal vector pointing in the +z direction.
1610 If the atom has z > 7, then a downward force would be applied of
1611 k * (atom.z - 7). The same plane with the normal vector pointing in
1612 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).
1613 """
1615 if isinstance(a2, int):
1616 self._type = 'two atoms'
1617 self.indices = [a1, a2]
1618 elif len(a2) == 3:
1619 self._type = 'point'
1620 self.index = a1
1621 self.origin = np.array(a2)
1622 elif len(a2) == 4:
1623 self._type = 'plane'
1624 self.index = a1
1625 self.plane = a2
1626 else:
1627 raise RuntimeError('Unknown type for a2')
1628 self.threshold = rt
1629 self.spring = k
1631 def get_removed_dof(self, atoms):
1632 return 0
1634 def todict(self):
1635 dct = {'name': 'Hookean'}
1636 dct['kwargs'] = {'rt': self.threshold,
1637 'k': self.spring}
1638 if self._type == 'two atoms':
1639 dct['kwargs']['a1'] = self.indices[0]
1640 dct['kwargs']['a2'] = self.indices[1]
1641 elif self._type == 'point':
1642 dct['kwargs']['a1'] = self.index
1643 dct['kwargs']['a2'] = self.origin
1644 elif self._type == 'plane':
1645 dct['kwargs']['a1'] = self.index
1646 dct['kwargs']['a2'] = self.plane
1647 else:
1648 raise NotImplementedError('Bad type: %s' % self._type)
1649 return dct
1651 def adjust_positions(self, atoms, newpositions):
1652 pass
1654 def adjust_momenta(self, atoms, momenta):
1655 pass
1657 def adjust_forces(self, atoms, forces):
1658 positions = atoms.positions
1659 if self._type == 'plane':
1660 A, B, C, D = self.plane
1661 x, y, z = positions[self.index]
1662 d = ((A * x + B * y + C * z + D) /
1663 np.sqrt(A**2 + B**2 + C**2))
1664 if d < 0:
1665 return
1666 magnitude = self.spring * d
1667 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))
1668 forces[self.index] += direction * magnitude
1669 return
1670 if self._type == 'two atoms':
1671 p1, p2 = positions[self.indices]
1672 elif self._type == 'point':
1673 p1 = positions[self.index]
1674 p2 = self.origin
1675 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1676 bondlength = np.linalg.norm(displace)
1677 if bondlength > self.threshold:
1678 magnitude = self.spring * (bondlength - self.threshold)
1679 direction = displace / np.linalg.norm(displace)
1680 if self._type == 'two atoms':
1681 forces[self.indices[0]] += direction * magnitude
1682 forces[self.indices[1]] -= direction * magnitude
1683 else:
1684 forces[self.index] += direction * magnitude
1686 def adjust_potential_energy(self, atoms):
1687 """Returns the difference to the potential energy due to an active
1688 constraint. (That is, the quantity returned is to be added to the
1689 potential energy.)"""
1690 positions = atoms.positions
1691 if self._type == 'plane':
1692 A, B, C, D = self.plane
1693 x, y, z = positions[self.index]
1694 d = ((A * x + B * y + C * z + D) /
1695 np.sqrt(A**2 + B**2 + C**2))
1696 if d > 0:
1697 return 0.5 * self.spring * d**2
1698 else:
1699 return 0.
1700 if self._type == 'two atoms':
1701 p1, p2 = positions[self.indices]
1702 elif self._type == 'point':
1703 p1 = positions[self.index]
1704 p2 = self.origin
1705 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)
1706 bondlength = np.linalg.norm(displace)
1707 if bondlength > self.threshold:
1708 return 0.5 * self.spring * (bondlength - self.threshold)**2
1709 else:
1710 return 0.
1712 def get_indices(self):
1713 if self._type == 'two atoms':
1714 return self.indices
1715 elif self._type == 'point':
1716 return self.index
1717 elif self._type == 'plane':
1718 return self.index
1720 def index_shuffle(self, atoms, ind):
1721 # See docstring of superclass
1722 if self._type == 'two atoms':
1723 newa = [-1, -1] # Signal error
1724 for new, old in slice2enlist(ind, len(atoms)):
1725 for i, a in enumerate(self.indices):
1726 if old == a:
1727 newa[i] = new
1728 if newa[0] == -1 or newa[1] == -1:
1729 raise IndexError('Constraint not part of slice')
1730 self.indices = newa
1731 elif (self._type == 'point') or (self._type == 'plane'):
1732 newa = -1 # Signal error
1733 for new, old in slice2enlist(ind, len(atoms)):
1734 if old == self.index:
1735 newa = new
1736 break
1737 if newa == -1:
1738 raise IndexError('Constraint not part of slice')
1739 self.index = newa
1741 def __repr__(self):
1742 if self._type == 'two atoms':
1743 return 'Hookean(%d, %d)' % tuple(self.indices)
1744 elif self._type == 'point':
1745 return 'Hookean(%d) to cartesian' % self.index
1746 else:
1747 return 'Hookean(%d) to plane' % self.index
1750class ExternalForce(FixConstraint):
1751 """Constraint object for pulling two atoms apart by an external force.
1753 You can combine this constraint for example with FixBondLength but make
1754 sure that *ExternalForce* comes first in the list if there are overlaps
1755 between atom1-2 and atom3-4:
1757 >>> con1 = ExternalForce(atom1, atom2, f_ext)
1758 >>> con2 = FixBondLength(atom3, atom4)
1759 >>> atoms.set_constraint([con1, con2])
1761 see ase/test/external_force.py"""
1763 def __init__(self, a1, a2, f_ext):
1764 self.indices = [a1, a2]
1765 self.external_force = f_ext
1767 def get_removed_dof(self, atoms):
1768 return 0
1770 def adjust_positions(self, atoms, new):
1771 pass
1773 def adjust_forces(self, atoms, forces):
1774 dist = np.subtract.reduce(atoms.positions[self.indices])
1775 force = self.external_force * dist / np.linalg.norm(dist)
1776 forces[self.indices] += (force, -force)
1778 def adjust_potential_energy(self, atoms):
1779 dist = np.subtract.reduce(atoms.positions[self.indices])
1780 return -np.linalg.norm(dist) * self.external_force
1782 def index_shuffle(self, atoms, ind):
1783 """Shuffle the indices of the two atoms in this constraint"""
1784 newa = [-1, -1] # Signal error
1785 for new, old in slice2enlist(ind, len(atoms)):
1786 for i, a in enumerate(self.indices):
1787 if old == a:
1788 newa[i] = new
1789 if newa[0] == -1 or newa[1] == -1:
1790 raise IndexError('Constraint not part of slice')
1791 self.indices = newa
1793 def __repr__(self):
1794 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],
1795 self.indices[1],
1796 self.external_force)
1798 def todict(self):
1799 return {'name': 'ExternalForce',
1800 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
1801 'f_ext': self.external_force}}
1804class MirrorForce(FixConstraint):
1805 """Constraint object for mirroring the force between two atoms.
1807 This class is designed to find a transition state with the help of a
1808 single optimization. It can be used if the transition state belongs to a
1809 bond breaking reaction. First the given bond length will be fixed until
1810 all other degrees of freedom are optimized, then the forces of the two
1811 atoms will be mirrored to find the transition state. The mirror plane is
1812 perpendicular to the connecting line of the atoms. Transition states in
1813 dependence of the force can be obtained by stretching the molecule and
1814 fixing its total length with *FixBondLength* or by using *ExternalForce*
1815 during the optimization with *MirrorForce*.
1817 Parameters
1818 ----------
1819 a1: int
1820 First atom index.
1821 a2: int
1822 Second atom index.
1823 max_dist: float
1824 Upper limit of the bond length interval where the transition state
1825 can be found.
1826 min_dist: float
1827 Lower limit of the bond length interval where the transition state
1828 can be found.
1829 fmax: float
1830 Maximum force used for the optimization.
1832 Notes
1833 -----
1834 You can combine this constraint for example with FixBondLength but make
1835 sure that *MirrorForce* comes first in the list if there are overlaps
1836 between atom1-2 and atom3-4:
1838 >>> con1 = MirrorForce(atom1, atom2)
1839 >>> con2 = FixBondLength(atom3, atom4)
1840 >>> atoms.set_constraint([con1, con2])
1842 """
1844 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):
1845 self.indices = [a1, a2]
1846 self.min_dist = min_dist
1847 self.max_dist = max_dist
1848 self.fmax = fmax
1850 def adjust_positions(self, atoms, new):
1851 pass
1853 def adjust_forces(self, atoms, forces):
1854 dist = np.subtract.reduce(atoms.positions[self.indices])
1855 d = np.linalg.norm(dist)
1856 if (d < self.min_dist) or (d > self.max_dist):
1857 # Stop structure optimization
1858 forces[:] *= 0
1859 return
1860 dist /= d
1861 df = np.subtract.reduce(forces[self.indices])
1862 f = df.dot(dist)
1863 con_saved = atoms.constraints
1864 try:
1865 con = [con for con in con_saved
1866 if not isinstance(con, MirrorForce)]
1867 atoms.set_constraint(con)
1868 forces_copy = atoms.get_forces()
1869 finally:
1870 atoms.set_constraint(con_saved)
1871 df1 = -1 / 2. * f * dist
1872 forces_copy[self.indices] += (df1, -df1)
1873 # Check if forces would be converged if the bond with mirrored forces
1874 # would also be fixed
1875 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
1876 factor = 1.
1877 else:
1878 factor = 0.
1879 df1 = -(1 + factor) / 2. * f * dist
1880 forces[self.indices] += (df1, -df1)
1882 def index_shuffle(self, atoms, ind):
1883 """Shuffle the indices of the two atoms in this constraint
1885 """
1886 newa = [-1, -1] # Signal error
1887 for new, old in slice2enlist(ind, len(atoms)):
1888 for i, a in enumerate(self.indices):
1889 if old == a:
1890 newa[i] = new
1891 if newa[0] == -1 or newa[1] == -1:
1892 raise IndexError('Constraint not part of slice')
1893 self.indices = newa
1895 def __repr__(self):
1896 return 'MirrorForce(%d, %d, %f, %f, %f)' % (
1897 self.indices[0], self.indices[1], self.max_dist, self.min_dist,
1898 self.fmax)
1900 def todict(self):
1901 return {'name': 'MirrorForce',
1902 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
1903 'max_dist': self.max_dist,
1904 'min_dist': self.min_dist, 'fmax': self.fmax}}
1907class MirrorTorque(FixConstraint):
1908 """Constraint object for mirroring the torque acting on a dihedral
1909 angle defined by four atoms.
1911 This class is designed to find a transition state with the help of a
1912 single optimization. It can be used if the transition state belongs to a
1913 cis-trans-isomerization with a change of dihedral angle. First the given
1914 dihedral angle will be fixed until all other degrees of freedom are
1915 optimized, then the torque acting on the dihedral angle will be mirrored
1916 to find the transition state. Transition states in
1917 dependence of the force can be obtained by stretching the molecule and
1918 fixing its total length with *FixBondLength* or by using *ExternalForce*
1919 during the optimization with *MirrorTorque*.
1921 This constraint can be used to find
1922 transition states of cis-trans-isomerization.
1924 a1 a4
1925 | |
1926 a2 __ a3
1928 Parameters
1929 ----------
1930 a1: int
1931 First atom index.
1932 a2: int
1933 Second atom index.
1934 a3: int
1935 Third atom index.
1936 a4: int
1937 Fourth atom index.
1938 max_angle: float
1939 Upper limit of the dihedral angle interval where the transition state
1940 can be found.
1941 min_angle: float
1942 Lower limit of the dihedral angle interval where the transition state
1943 can be found.
1944 fmax: float
1945 Maximum force used for the optimization.
1947 Notes
1948 -----
1949 You can combine this constraint for example with FixBondLength but make
1950 sure that *MirrorTorque* comes first in the list if there are overlaps
1951 between atom1-4 and atom5-6:
1953 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)
1954 >>> con2 = FixBondLength(atom5, atom6)
1955 >>> atoms.set_constraint([con1, con2])
1957 """
1959 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,
1960 fmax=0.1):
1961 self.indices = [a1, a2, a3, a4]
1962 self.min_angle = min_angle
1963 self.max_angle = max_angle
1964 self.fmax = fmax
1966 def adjust_positions(self, atoms, new):
1967 pass
1969 def adjust_forces(self, atoms, forces):
1970 angle = atoms.get_dihedral(self.indices[0], self.indices[1],
1971 self.indices[2], self.indices[3])
1972 angle *= np.pi / 180.
1973 if (angle < self.min_angle) or (angle > self.max_angle):
1974 # Stop structure optimization
1975 forces[:] *= 0
1976 return
1977 p = atoms.positions[self.indices]
1978 f = forces[self.indices]
1980 f0 = (f[1] + f[2]) / 2.
1981 ff = f - f0
1982 p0 = (p[2] + p[1]) / 2.
1983 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)
1984 fff = ff - np.cross(m0, p - p0)
1985 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \
1986 (p[1] - p0).dot(p[1] - p0)
1987 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \
1988 (p[2] - p0).dot(p[2] - p0)
1989 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \
1990 np.linalg.norm(p[1] - p0)
1991 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \
1992 np.linalg.norm(p[2] - p0)
1993 omegap = omegap1 + omegap2
1994 con_saved = atoms.constraints
1995 try:
1996 con = [con for con in con_saved
1997 if not isinstance(con, MirrorTorque)]
1998 atoms.set_constraint(con)
1999 forces_copy = atoms.get_forces()
2000 finally:
2001 atoms.set_constraint(con_saved)
2002 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \
2003 np.linalg.norm(p[1] - p0)
2004 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \
2005 np.linalg.norm(p[2] - p0)
2006 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2007 # Check if forces would be converged if the dihedral angle with
2008 # mirrored torque would also be fixed
2009 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:
2010 factor = 1.
2011 else:
2012 factor = 0.
2013 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \
2014 np.linalg.norm(p[1] - p0)
2015 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \
2016 np.linalg.norm(p[2] - p0)
2017 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)
2019 def index_shuffle(self, atoms, ind):
2020 # See docstring of superclass
2021 indices = []
2022 for new, old in slice2enlist(ind, len(atoms)):
2023 if old in self.indices:
2024 indices.append(new)
2025 if len(indices) == 0:
2026 raise IndexError('All indices in MirrorTorque not part of slice')
2027 self.indices = np.asarray(indices, int)
2029 def __repr__(self):
2030 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (
2031 self.indices[0], self.indices[1], self.indices[2],
2032 self.indices[3], self.max_angle, self.min_angle, self.fmax)
2034 def todict(self):
2035 return {'name': 'MirrorTorque',
2036 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],
2037 'a3': self.indices[2], 'a4': self.indices[3],
2038 'max_angle': self.max_angle,
2039 'min_angle': self.min_angle, 'fmax': self.fmax}}
2042class Filter:
2043 """Subset filter class."""
2045 def __init__(self, atoms, indices=None, mask=None):
2046 """Filter atoms.
2048 This filter can be used to hide degrees of freedom in an Atoms
2049 object.
2051 Parameters
2052 ----------
2053 indices : list of int
2054 Indices for those atoms that should remain visible.
2055 mask : list of bool
2056 One boolean per atom indicating if the atom should remain
2057 visible or not.
2059 If a Trajectory tries to save this object, it will instead
2060 save the underlying Atoms object. To prevent this, override
2061 the iterimages method.
2062 """
2064 self.atoms = atoms
2065 self.constraints = []
2066 # Make self.info a reference to the underlying atoms' info dictionary.
2067 self.info = self.atoms.info
2069 if indices is None and mask is None:
2070 raise ValueError('Use "indices" or "mask".')
2071 if indices is not None and mask is not None:
2072 raise ValueError('Use only one of "indices" and "mask".')
2074 if mask is not None:
2075 self.index = np.asarray(mask, bool)
2076 self.n = self.index.sum()
2077 else:
2078 self.index = np.asarray(indices, int)
2079 self.n = len(self.index)
2081 def iterimages(self):
2082 # Present the real atoms object to Trajectory and friends
2083 return self.atoms.iterimages()
2085 def get_cell(self):
2086 """Returns the computational cell.
2088 The computational cell is the same as for the original system.
2089 """
2090 return self.atoms.get_cell()
2092 def get_pbc(self):
2093 """Returns the periodic boundary conditions.
2095 The boundary conditions are the same as for the original system.
2096 """
2097 return self.atoms.get_pbc()
2099 def get_positions(self):
2100 'Return the positions of the visible atoms.'
2101 return self.atoms.get_positions()[self.index]
2103 def set_positions(self, positions, **kwargs):
2104 'Set the positions of the visible atoms.'
2105 pos = self.atoms.get_positions()
2106 pos[self.index] = positions
2107 self.atoms.set_positions(pos, **kwargs)
2109 positions = property(get_positions, set_positions,
2110 doc='Positions of the atoms')
2112 def get_momenta(self):
2113 'Return the momenta of the visible atoms.'
2114 return self.atoms.get_momenta()[self.index]
2116 def set_momenta(self, momenta, **kwargs):
2117 'Set the momenta of the visible atoms.'
2118 mom = self.atoms.get_momenta()
2119 mom[self.index] = momenta
2120 self.atoms.set_momenta(mom, **kwargs)
2122 def get_atomic_numbers(self):
2123 'Return the atomic numbers of the visible atoms.'
2124 return self.atoms.get_atomic_numbers()[self.index]
2126 def set_atomic_numbers(self, atomic_numbers):
2127 'Set the atomic numbers of the visible atoms.'
2128 z = self.atoms.get_atomic_numbers()
2129 z[self.index] = atomic_numbers
2130 self.atoms.set_atomic_numbers(z)
2132 def get_tags(self):
2133 'Return the tags of the visible atoms.'
2134 return self.atoms.get_tags()[self.index]
2136 def set_tags(self, tags):
2137 'Set the tags of the visible atoms.'
2138 tg = self.atoms.get_tags()
2139 tg[self.index] = tags
2140 self.atoms.set_tags(tg)
2142 def get_forces(self, *args, **kwargs):
2143 return self.atoms.get_forces(*args, **kwargs)[self.index]
2145 def get_stress(self, *args, **kwargs):
2146 return self.atoms.get_stress(*args, **kwargs)
2148 def get_stresses(self, *args, **kwargs):
2149 return self.atoms.get_stresses(*args, **kwargs)[self.index]
2151 def get_masses(self):
2152 return self.atoms.get_masses()[self.index]
2154 def get_potential_energy(self, **kwargs):
2155 """Calculate potential energy.
2157 Returns the potential energy of the full system.
2158 """
2159 return self.atoms.get_potential_energy(**kwargs)
2161 def get_chemical_symbols(self):
2162 return self.atoms.get_chemical_symbols()
2164 def get_initial_magnetic_moments(self):
2165 return self.atoms.get_initial_magnetic_moments()
2167 def get_calculator(self):
2168 """Returns the calculator.
2170 WARNING: The calculator is unaware of this filter, and sees a
2171 different number of atoms.
2172 """
2173 return self.atoms.calc
2175 @property
2176 def calc(self):
2177 return self.atoms.calc
2179 def get_celldisp(self):
2180 return self.atoms.get_celldisp()
2182 def has(self, name):
2183 'Check for existence of array.'
2184 return self.atoms.has(name)
2186 def __len__(self):
2187 'Return the number of movable atoms.'
2188 return self.n
2190 def __getitem__(self, i):
2191 'Return an atom.'
2192 return self.atoms[self.index[i]]
2195class StrainFilter(Filter):
2196 """Modify the supercell while keeping the scaled positions fixed.
2198 Presents the strain of the supercell as the generalized positions,
2199 and the global stress tensor (times the volume) as the generalized
2200 force.
2202 This filter can be used to relax the unit cell until the stress is
2203 zero. If MDMin is used for this, the timestep (dt) to be used
2204 depends on the system size. 0.01/x where x is a typical dimension
2205 seems like a good choice.
2207 The stress and strain are presented as 6-vectors, the order of the
2208 components follow the standard engingeering practice: xx, yy, zz,
2209 yz, xz, xy.
2211 """
2213 def __init__(self, atoms, mask=None, include_ideal_gas=False):
2214 """Create a filter applying a homogeneous strain to a list of atoms.
2216 The first argument, atoms, is the atoms object.
2218 The optional second argument, mask, is a list of six booleans,
2219 indicating which of the six independent components of the
2220 strain that are allowed to become non-zero. It defaults to
2221 [1,1,1,1,1,1].
2223 """
2225 self.strain = np.zeros(6)
2226 self.include_ideal_gas = include_ideal_gas
2228 if mask is None:
2229 mask = np.ones(6)
2230 else:
2231 mask = np.array(mask)
2233 Filter.__init__(self, atoms, mask=mask)
2234 self.mask = mask
2235 self.origcell = atoms.get_cell()
2237 def get_positions(self):
2238 return self.strain.reshape((2, 3)).copy()
2240 def set_positions(self, new):
2241 new = new.ravel() * self.mask
2242 eps = np.array([[1.0 + new[0], 0.5 * new[5], 0.5 * new[4]],
2243 [0.5 * new[5], 1.0 + new[1], 0.5 * new[3]],
2244 [0.5 * new[4], 0.5 * new[3], 1.0 + new[2]]])
2246 self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True)
2247 self.strain[:] = new
2249 def get_forces(self, **kwargs):
2250 stress = self.atoms.get_stress(include_ideal_gas=self.include_ideal_gas)
2251 return -self.atoms.get_volume() * (stress * self.mask).reshape((2, 3))
2253 def has(self, x):
2254 return self.atoms.has(x)
2256 def __len__(self):
2257 return 2
2260class UnitCellFilter(Filter):
2261 """Modify the supercell and the atom positions. """
2262 def __init__(self, atoms, mask=None,
2263 cell_factor=None,
2264 hydrostatic_strain=False,
2265 constant_volume=False,
2266 scalar_pressure=0.0):
2267 """Create a filter that returns the atomic forces and unit cell
2268 stresses together, so they can simultaneously be minimized.
2270 The first argument, atoms, is the atoms object. The optional second
2271 argument, mask, is a list of booleans, indicating which of the six
2272 independent components of the strain are relaxed.
2274 - True = relax to zero
2275 - False = fixed, ignore this component
2277 Degrees of freedom are the positions in the original undeformed cell,
2278 plus the deformation tensor (extra 3 "atoms"). This gives forces
2279 consistent with numerical derivatives of the potential energy
2280 with respect to the cell degreees of freedom.
2282 For full details see:
2283 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
2284 Phys. Rev. B 59, 235 (1999)
2286 You can still use constraints on the atoms, e.g. FixAtoms, to control
2287 the relaxation of the atoms.
2289 >>> # this should be equivalent to the StrainFilter
2290 >>> atoms = Atoms(...)
2291 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
2292 >>> ucf = UnitCellFilter(atoms)
2294 You should not attach this UnitCellFilter object to a
2295 trajectory. Instead, create a trajectory for the atoms, and
2296 attach it to an optimizer like this:
2298 >>> atoms = Atoms(...)
2299 >>> ucf = UnitCellFilter(atoms)
2300 >>> qn = QuasiNewton(ucf)
2301 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
2302 >>> qn.attach(traj)
2303 >>> qn.run(fmax=0.05)
2305 Helpful conversion table:
2307 - 0.05 eV/A^3 = 8 GPA
2308 - 0.003 eV/A^3 = 0.48 GPa
2309 - 0.0006 eV/A^3 = 0.096 GPa
2310 - 0.0003 eV/A^3 = 0.048 GPa
2311 - 0.0001 eV/A^3 = 0.02 GPa
2313 Additional optional arguments:
2315 cell_factor: float (default float(len(atoms)))
2316 Factor by which deformation gradient is multiplied to put
2317 it on the same scale as the positions when assembling
2318 the combined position/cell vector. The stress contribution to
2319 the forces is scaled down by the same factor. This can be thought
2320 of as a very simple preconditioners. Default is number of atoms
2321 which gives approximately the correct scaling.
2323 hydrostatic_strain: bool (default False)
2324 Constrain the cell by only allowing hydrostatic deformation.
2325 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
2327 constant_volume: bool (default False)
2328 Project out the diagonal elements of the virial tensor to allow
2329 relaxations at constant volume, e.g. for mapping out an
2330 energy-volume curve. Note: this only approximately conserves
2331 the volume and breaks energy/force consistency so can only be
2332 used with optimizers that do require do a line minimisation
2333 (e.g. FIRE).
2335 scalar_pressure: float (default 0.0)
2336 Applied pressure to use for enthalpy pV term. As above, this
2337 breaks energy/force consistency.
2338 """
2340 Filter.__init__(self, atoms, indices=range(len(atoms)))
2341 self.atoms = atoms
2342 self.orig_cell = atoms.get_cell()
2343 self.stress = None
2345 if mask is None:
2346 mask = np.ones(6)
2347 mask = np.asarray(mask)
2348 if mask.shape == (6,):
2349 self.mask = voigt_6_to_full_3x3_stress(mask)
2350 elif mask.shape == (3, 3):
2351 self.mask = mask
2352 else:
2353 raise ValueError('shape of mask should be (3,3) or (6,)')
2355 if cell_factor is None:
2356 cell_factor = float(len(atoms))
2357 self.hydrostatic_strain = hydrostatic_strain
2358 self.constant_volume = constant_volume
2359 self.scalar_pressure = scalar_pressure
2360 self.cell_factor = cell_factor
2361 self.copy = self.atoms.copy
2362 self.arrays = self.atoms.arrays
2364 def deform_grad(self):
2365 return np.linalg.solve(self.orig_cell, self.atoms.cell).T
2367 def get_positions(self):
2368 """
2369 this returns an array with shape (natoms + 3,3).
2371 the first natoms rows are the positions of the atoms, the last
2372 three rows are the deformation tensor associated with the unit cell,
2373 scaled by self.cell_factor.
2374 """
2376 cur_deform_grad = self.deform_grad()
2377 natoms = len(self.atoms)
2378 pos = np.zeros((natoms + 3, 3))
2379 # UnitCellFilter's positions are the self.atoms.positions but without
2380 # the applied deformation gradient
2381 pos[:natoms] = np.linalg.solve(cur_deform_grad,
2382 self.atoms.positions.T).T
2383 # UnitCellFilter's cell DOFs are the deformation gradient times a
2384 # scaling factor
2385 pos[natoms:] = self.cell_factor * cur_deform_grad
2386 return pos
2388 def set_positions(self, new, **kwargs):
2389 """
2390 new is an array with shape (natoms+3,3).
2392 the first natoms rows are the positions of the atoms, the last
2393 three rows are the deformation tensor used to change the cell shape.
2395 the new cell is first set from original cell transformed by the new
2396 deformation gradient, then the positions are set with respect to the
2397 current cell by transforming them with the same deformation gradient
2398 """
2400 natoms = len(self.atoms)
2401 new_atom_positions = new[:natoms]
2402 new_deform_grad = new[natoms:] / self.cell_factor
2403 # Set the new cell from the original cell and the new
2404 # deformation gradient. Both current and final structures should
2405 # preserve symmetry, so if set_cell() calls FixSymmetry.adjust_cell(),
2406 # it should be OK
2407 self.atoms.set_cell(self.orig_cell @ new_deform_grad.T,
2408 scale_atoms=True)
2409 # Set the positions from the ones passed in (which are without the
2410 # deformation gradient applied) and the new deformation gradient.
2411 # This should also preserve symmetry, so if set_positions() calls
2412 # FixSymmetyr.adjust_positions(), it should be OK
2413 self.atoms.set_positions(new_atom_positions @ new_deform_grad.T,
2414 **kwargs)
2416 def get_potential_energy(self, force_consistent=True):
2417 """
2418 returns potential energy including enthalpy PV term.
2419 """
2420 atoms_energy = self.atoms.get_potential_energy(
2421 force_consistent=force_consistent)
2422 return atoms_energy + self.scalar_pressure * self.atoms.get_volume()
2424 def get_forces(self, **kwargs):
2425 """
2426 returns an array with shape (natoms+3,3) of the atomic forces
2427 and unit cell stresses.
2429 the first natoms rows are the forces on the atoms, the last
2430 three rows are the forces on the unit cell, which are
2431 computed from the stress tensor.
2432 """
2434 stress = self.atoms.get_stress(**kwargs)
2435 atoms_forces = self.atoms.get_forces(**kwargs)
2437 volume = self.atoms.get_volume()
2438 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
2439 np.diag([self.scalar_pressure] * 3))
2440 cur_deform_grad = self.deform_grad()
2441 atoms_forces = atoms_forces @ cur_deform_grad
2442 virial = np.linalg.solve(cur_deform_grad, virial.T).T
2444 if self.hydrostatic_strain:
2445 vtr = virial.trace()
2446 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
2448 # Zero out components corresponding to fixed lattice elements
2449 if (self.mask != 1.0).any():
2450 virial *= self.mask
2452 if self.constant_volume:
2453 vtr = virial.trace()
2454 np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0)
2456 natoms = len(self.atoms)
2457 forces = np.zeros((natoms + 3, 3))
2458 forces[:natoms] = atoms_forces
2459 forces[natoms:] = virial / self.cell_factor
2461 self.stress = -full_3x3_to_voigt_6_stress(virial) / volume
2462 return forces
2464 def get_stress(self):
2465 raise PropertyNotImplementedError
2467 def has(self, x):
2468 return self.atoms.has(x)
2470 def __len__(self):
2471 return (len(self.atoms) + 3)
2474class ExpCellFilter(UnitCellFilter):
2475 """Modify the supercell and the atom positions."""
2476 def __init__(self, atoms, mask=None,
2477 cell_factor=None,
2478 hydrostatic_strain=False,
2479 constant_volume=False,
2480 scalar_pressure=0.0):
2481 r"""Create a filter that returns the atomic forces and unit cell
2482 stresses together, so they can simultaneously be minimized.
2484 The first argument, atoms, is the atoms object. The optional second
2485 argument, mask, is a list of booleans, indicating which of the six
2486 independent components of the strain are relaxed.
2488 - True = relax to zero
2489 - False = fixed, ignore this component
2491 Degrees of freedom are the positions in the original undeformed cell,
2492 plus the log of the deformation tensor (extra 3 "atoms"). This gives
2493 forces consistent with numerical derivatives of the potential energy
2494 with respect to the cell degrees of freedom.
2496 For full details see:
2497 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
2498 Phys. Rev. B 59, 235 (1999)
2500 You can still use constraints on the atoms, e.g. FixAtoms, to control
2501 the relaxation of the atoms.
2503 >>> # this should be equivalent to the StrainFilter
2504 >>> atoms = Atoms(...)
2505 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
2506 >>> ecf = ExpCellFilter(atoms)
2508 You should not attach this ExpCellFilter object to a
2509 trajectory. Instead, create a trajectory for the atoms, and
2510 attach it to an optimizer like this:
2512 >>> atoms = Atoms(...)
2513 >>> ecf = ExpCellFilter(atoms)
2514 >>> qn = QuasiNewton(ecf)
2515 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
2516 >>> qn.attach(traj)
2517 >>> qn.run(fmax=0.05)
2519 Helpful conversion table:
2521 - 0.05 eV/A^3 = 8 GPA
2522 - 0.003 eV/A^3 = 0.48 GPa
2523 - 0.0006 eV/A^3 = 0.096 GPa
2524 - 0.0003 eV/A^3 = 0.048 GPa
2525 - 0.0001 eV/A^3 = 0.02 GPa
2527 Additional optional arguments:
2529 cell_factor: (DEPRECATED)
2530 Retained for backwards compatibility, but no longer used.
2532 hydrostatic_strain: bool (default False)
2533 Constrain the cell by only allowing hydrostatic deformation.
2534 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
2536 constant_volume: bool (default False)
2537 Project out the diagonal elements of the virial tensor to allow
2538 relaxations at constant volume, e.g. for mapping out an
2539 energy-volume curve.
2541 scalar_pressure: float (default 0.0)
2542 Applied pressure to use for enthalpy pV term. As above, this
2543 breaks energy/force consistency.
2545 Implementation details:
2547 The implementation is based on that of Christoph Ortner in JuLIP.jl:
2548 https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244
2550 We decompose the deformation gradient as
2552 F = exp(U) F0
2553 x = F * F0^{-1} z = exp(U) z
2555 If we write the energy as a function of U we can transform the
2556 stress associated with a perturbation V into a derivative using a
2557 linear map V -> L(U, V).
2559 \phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) .
2560 ( L(U, V) exp(-U) exp(U) z )
2562 where
2564 \nabla E(U) : V = [S exp(-U)'] : L(U,V)
2565 = L'(U, S exp(-U)') : V
2566 = L(U', S exp(-U)') : V
2567 = L(U, S exp(-U)) : V (provided U = U')
2569 where the : operator represents double contraction,
2570 i.e. A:B = trace(A'B), and
2572 F = deformation tensor - 3x3 matrix
2573 F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here
2574 U = cell degrees of freedom used here - 3x3 matrix
2575 V = perturbation to cell DoFs - 3x3 matrix
2576 v = perturbation to position DoFs
2577 x = atomic positions in deformed cell
2578 z = atomic positions in original cell
2579 \phi = potential energy
2580 S = stress tensor [3x3 matrix]
2581 L(U, V) = directional derivative of exp at U in direction V, i.e
2582 d/dt exp(U + t V)|_{t=0} = L(U, V)
2584 This means we can write
2586 d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V
2588 and therefore the contribution to the gradient of the energy is
2590 \nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij
2591 """
2593 Filter.__init__(self, atoms, indices=range(len(atoms)))
2594 UnitCellFilter.__init__(self, atoms, mask, cell_factor,
2595 hydrostatic_strain,
2596 constant_volume, scalar_pressure)
2597 if cell_factor is not None:
2598 warn("cell_factor is deprecated")
2599 self.cell_factor = 1.0
2601 def get_positions(self):
2602 pos = UnitCellFilter.get_positions(self)
2603 natoms = len(self.atoms)
2604 pos[natoms:] = logm(self.deform_grad())
2605 return pos
2607 def set_positions(self, new, **kwargs):
2608 natoms = len(self.atoms)
2609 new2 = new.copy()
2610 new2[natoms:] = expm(new[natoms:])
2611 UnitCellFilter.set_positions(self, new2, **kwargs)
2613 def get_forces(self, **kwargs):
2614 forces = UnitCellFilter.get_forces(self, **kwargs)
2616 # forces on atoms are same as UnitCellFilter, we just
2617 # need to modify the stress contribution
2618 stress = self.atoms.get_stress(**kwargs)
2619 volume = self.atoms.get_volume()
2620 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
2621 np.diag([self.scalar_pressure] * 3))
2623 cur_deform_grad = self.deform_grad()
2624 cur_deform_grad_log = logm(cur_deform_grad)
2626 if self.hydrostatic_strain:
2627 vtr = virial.trace()
2628 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
2630 # Zero out components corresponding to fixed lattice elements
2631 if (self.mask != 1.0).any():
2632 virial *= self.mask
2634 deform_grad_log_force_naive = virial.copy()
2635 Y = np.zeros((6, 6))
2636 Y[0:3, 0:3] = cur_deform_grad_log
2637 Y[3:6, 3:6] = cur_deform_grad_log
2638 Y[0:3, 3:6] = - virial @ expm(-cur_deform_grad_log)
2639 deform_grad_log_force = -expm(Y)[0:3, 3:6]
2640 for (i1, i2) in [(0, 1), (0, 2), (1, 2)]:
2641 ff = 0.5 * (deform_grad_log_force[i1, i2] +
2642 deform_grad_log_force[i2, i1])
2643 deform_grad_log_force[i1, i2] = ff
2644 deform_grad_log_force[i2, i1] = ff
2646 # check for reasonable alignment between naive and
2647 # exact search directions
2648 all_are_equal = np.all(np.isclose(deform_grad_log_force,
2649 deform_grad_log_force_naive))
2650 if all_are_equal or \
2651 (np.sum(deform_grad_log_force * deform_grad_log_force_naive) /
2652 np.sqrt(np.sum(deform_grad_log_force**2) *
2653 np.sum(deform_grad_log_force_naive**2)) > 0.8):
2654 deform_grad_log_force = deform_grad_log_force_naive
2656 # Cauchy stress used for convergence testing
2657 convergence_crit_stress = -(virial / volume)
2658 if self.constant_volume:
2659 # apply constraint to force
2660 dglf_trace = deform_grad_log_force.trace()
2661 np.fill_diagonal(deform_grad_log_force,
2662 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
2663 # apply constraint to Cauchy stress used for convergence testing
2664 ccs_trace = convergence_crit_stress.trace()
2665 np.fill_diagonal(convergence_crit_stress,
2666 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
2668 # pack gradients into vector
2669 natoms = len(self.atoms)
2670 forces[natoms:] = deform_grad_log_force
2671 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
2672 return forces