1from math import sqrt 2from warnings import warn 3 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) 13 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"] 21 22 23def dict2constraint(dct): 24 if dct['name'] not in __all__: 25 raise ValueError 26 return globals()[dct['name']](**dct['kwargs']) 27 28 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) 34 35 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)) 49 50 51class FixConstraint: 52 """Base class for classes that fix one or more atoms in some way.""" 53 54 def index_shuffle(self, atoms, ind): 55 """Change the indices. 56 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. 60 61 ind -- List or tuple of indices. 62 63 """ 64 raise NotImplementedError 65 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) 75 76 def adjust_momenta(self, atoms, momenta): 77 """Adjusts momenta in identical manner to forces.""" 78 self.adjust_forces(atoms, momenta) 79 80 def copy(self): 81 return dict2constraint(self.todict().copy()) 82 83 84class FixConstraintSingle(FixConstraint): 85 """Base class for classes that fix a single atom.""" 86 87 def __init__(self, a): 88 self.a = a 89 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 102 103 def get_indices(self): 104 return [self.a] 105 106 107class FixAtoms(FixConstraint): 108 """Constraint object for fixing some chosen atoms.""" 109 110 def __init__(self, indices=None, mask=None): 111 """Constrain chosen atoms. 112 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. 120 121 Examples 122 -------- 123 Fix all Copper atoms: 124 125 >>> mask = [s == 'Cu' for s in atoms.get_chemical_symbols()] 126 >>> c = FixAtoms(mask=mask) 127 >>> atoms.set_constraint(c) 128 129 Fix all atoms with z-coordinate less than 1.0 Angstrom: 130 131 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0) 132 >>> atoms.set_constraint(c) 133 """ 134 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".') 139 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) 151 152 if self.index.ndim != 1: 153 raise ValueError('Wrong argument to FixAtoms class!') 154 155 def get_removed_dof(self, atoms): 156 return 3 * len(self.index) 157 158 def adjust_positions(self, atoms, new): 159 new[self.index] = atoms.positions[self.index] 160 161 def adjust_forces(self, atoms, forces): 162 forces[self.index] = 0.0 163 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) 173 174 def get_indices(self): 175 return self.index 176 177 def __repr__(self): 178 return 'FixAtoms(indices=%s)' % ints2string(self.index) 179 180 def todict(self): 181 return {'name': 'FixAtoms', 182 'kwargs': {'indices': self.index.tolist()}} 183 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 199 200 def delete_atoms(self, indices, natoms): 201 """Removes atom number ind from the index array, if present. 202 203 Required for removing atoms with existing FixAtoms constraints. 204 """ 205 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 214 215 216class FixCom(FixConstraint): 217 """Constraint class for fixing the center of mass. 218 219 References 220 221 https://pubs.acs.org/doi/abs/10.1021/jp9722824 222 223 """ 224 225 def get_removed_dof(self, atoms): 226 return 3 227 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 234 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 240 241 def todict(self): 242 return {'name': 'FixCom', 243 'kwargs': {}} 244 245 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] + ', ...]' 251 252 253class FixBondLengths(FixConstraint): 254 maxiter = 500 255 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 263 264 def get_removed_dof(self, atoms): 265 return len(self.pairs) 266 267 def adjust_positions(self, atoms, new): 268 old = atoms.positions 269 masses = atoms.get_masses() 270 271 if self.bondlengths is None: 272 self.bondlengths = self.initialize_bond_lengths(atoms) 273 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') 293 294 def adjust_momenta(self, atoms, p): 295 old = atoms.positions 296 masses = atoms.get_masses() 297 298 if self.bondlengths is None: 299 self.bondlengths = self.initialize_bond_lengths(atoms) 300 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') 320 321 def adjust_forces(self, atoms, forces): 322 self.constraint_forces = -forces 323 self.adjust_momenta(atoms, forces) 324 self.constraint_forces += forces 325 326 def initialize_bond_lengths(self, atoms): 327 bondlengths = np.zeros(len(self.pairs)) 328 329 for i, ab in enumerate(self.pairs): 330 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True) 331 332 return bondlengths 333 334 def get_indices(self): 335 return np.unique(self.pairs.ravel()) 336 337 def todict(self): 338 return {'name': 'FixBondLengths', 339 'kwargs': {'pairs': self.pairs.tolist(), 340 'tolerance': self.tolerance}} 341 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') 353 354 355def FixBondLength(a1, a2): 356 """Fix distance between atoms with indices a1 and a2.""" 357 return FixBondLengths([(a1, a2)]) 358 359 360class FixLinearTriatomic(FixConstraint): 361 """Holonomic constraints for rigid linear triatomic molecules.""" 362 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: 368 369 n--o--m 370 371 Parameters: 372 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). 376 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. 383 384 References: 385 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 393 394 def get_removed_dof(self, atoms): 395 return 4 * len(self.triples) 396 397 @property 398 def n_ind(self): 399 return self.triples[:, 0] 400 401 @property 402 def m_ind(self): 403 return self.triples[:, 2] 404 405 @property 406 def o_ind(self): 407 return self.triples[:, 1] 408 409 def initialize(self, atoms): 410 masses = atoms.get_masses() 411 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses) 412 413 self.bondlengths = self.initialize_bond_lengths(atoms) 414 self.bondlengths_nm = self.bondlengths.sum(axis=1) 415 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] 427 428 self.C1 = C1 429 self.C2 = C2 430 self.C3 = C3 431 self.C4 = C4 432 433 def adjust_positions(self, atoms, new): 434 old = atoms.positions 435 new_n, new_m, new_o = self.get_slices(new) 436 437 if self.bondlengths is None: 438 self.initialize(atoms) 439 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) 458 459 self.set_slices(new_n, new_m, new_o, new) 460 461 def adjust_momenta(self, atoms, p): 462 old = atoms.positions 463 p_n, p_m, p_o = self.get_slices(p) 464 465 if self.bondlengths is None: 466 self.initialize(atoms) 467 468 mass_nn = self.mass_n[:, None] 469 mass_mm = self.mass_m[:, None] 470 mass_oo = self.mass_o[:, None] 471 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)) 481 482 self.set_slices(p_n, p_m, p_o, p) 483 484 def adjust_forces(self, atoms, forces): 485 486 if self.bondlengths is None: 487 self.initialize(atoms) 488 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] 494 495 self.constraint_forces = -forces 496 old = atoms.positions 497 498 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces) 499 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 507 508 self.constraint_forces += forces 509 510 def redistribute_forces_optimization(self, forces): 511 """Redistribute forces within a triple when performing structure 512 optimizations. 513 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] 521 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) 528 529 return fr_n, fr_m, fr_o 530 531 def redistribute_forces_md(self, atoms, forces, rand=False): 532 """Redistribute forces within a triple when performing molecular 533 dynamics. 534 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 559 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)) 563 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)) 567 568 self.set_slices(fr_n, fr_m, 0.0, forces) 569 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] 574 575 return a_n, a_m, a_o 576 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 581 582 def initialize_bond_lengths(self, atoms): 583 bondlengths = np.zeros((len(self.triples), 2)) 584 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) 590 591 return bondlengths 592 593 def get_indices(self): 594 return np.unique(self.triples.ravel()) 595 596 def todict(self): 597 return {'name': 'FixLinearTriatomic', 598 'kwargs': {'triples': self.triples.tolist()}} 599 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') 611 612 613class FixedMode(FixConstraint): 614 """Constrain atoms to move along directions orthogonal to 615 a given mode only.""" 616 617 def __init__(self, mode): 618 self.mode = (np.asarray(mode) / np.sqrt((mode**2).sum())).reshape(-1) 619 620 def get_removed_dof(self, atoms): 621 return len(atoms) 622 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) 628 629 def adjust_forces(self, atoms, forces): 630 forces = forces.ravel() 631 forces -= self.mode * np.dot(forces, self.mode) 632 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() 641 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 [] 648 649 def todict(self): 650 return {'name': 'FixedMode', 651 'kwargs': {'mode': self.mode.tolist()}} 652 653 def __repr__(self): 654 return 'FixedMode(%s)' % self.mode.tolist() 655 656 657class FixedPlane(FixConstraintSingle): 658 """Constrain an atom index *a* to move in a given plane only. 659 660 The plane is defined by its normal vector *direction*.""" 661 662 def __init__(self, a, direction): 663 self.a = a 664 self.dir = np.asarray(direction) / sqrt(np.dot(direction, direction)) 665 666 def get_removed_dof(self, atoms): 667 return 1 668 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) 672 673 def adjust_forces(self, atoms, forces): 674 forces[self.a] -= self.dir * np.dot(forces[self.a], self.dir) 675 676 def todict(self): 677 return {'name': 'FixedPlane', 678 'kwargs': {'a': self.a, 'direction': self.dir.tolist()}} 679 680 def __repr__(self): 681 return 'FixedPlane(%d, %s)' % (self.a, self.dir.tolist()) 682 683 684class FixedLine(FixConstraintSingle): 685 """Constrain an atom index *a* to move on a given line only. 686 687 The line is defined by its vector *direction*.""" 688 689 def __init__(self, a, direction): 690 self.a = a 691 self.dir = np.asarray(direction) / sqrt(np.dot(direction, direction)) 692 693 def get_removed_dof(self, atoms): 694 return 2 695 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 700 701 def adjust_forces(self, atoms, forces): 702 forces[self.a] = self.dir * np.dot(forces[self.a], self.dir) 703 704 def __repr__(self): 705 return 'FixedLine(%d, %s)' % (self.a, self.dir.tolist()) 706 707 def todict(self): 708 return {'name': 'FixedLine', 709 'kwargs': {'a': self.a, 'direction': self.dir.tolist()}} 710 711 712class FixCartesian(FixConstraintSingle): 713 'Fix an atom index *a* in the directions of the cartesian coordinates.' 714 715 def __init__(self, a, mask=(1, 1, 1)): 716 self.a = a 717 self.mask = ~np.asarray(mask, bool) 718 719 def get_removed_dof(self, atoms): 720 return 3 - self.mask.sum() 721 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 726 727 def adjust_forces(self, atoms, forces): 728 forces[self.a] *= self.mask 729 730 def __repr__(self): 731 return 'FixCartesian(a={0}, mask={1})'.format(self.a, 732 list(~self.mask)) 733 734 def todict(self): 735 return {'name': 'FixCartesian', 736 'kwargs': {'a': self.a, 'mask': ~self.mask.tolist()}} 737 738 739class FixScaled(FixConstraintSingle): 740 'Fix an atom index *a* in the directions of the unit vectors.' 741 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) 746 747 def get_removed_dof(self, atoms): 748 return self.mask.sum() 749 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] 757 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] 764 765 def todict(self): 766 return {'name': 'FixScaled', 767 'kwargs': {'a': self.a, 768 'cell': self.cell.tolist(), 769 'mask': self.mask.tolist()}} 770 771 def __repr__(self): 772 return 'FixScaled(%s, %d, %s)' % (repr(self.cell), 773 self.a, 774 repr(self.mask)) 775 776 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. 781 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): 790 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 807 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 816 817 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals) 818 + len(self.bondcombos) + len(self.anglecombos) 819 + len(self.dihedralcombos)) 820 821 # Initialize these at run-time: 822 self.constraints = [] 823 self.initialized = False 824 825 def get_removed_dof(self, atoms): 826 return self.n 827 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 848 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 863 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 882 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') 899 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)) 908 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}} 919 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.') 932 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 940 941 list_constraints = [r.ravel() for r in list2_constraints] 942 943 tx[:, 0] = 1.0 944 ty[:, 1] = 1.0 945 tz[:, 2] = 1.0 946 ff = forces.ravel() 947 948 # Calculate the center of mass 949 center = positions.sum(axis=0) / N 950 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] 957 958 # Normalizing transl., rotat. constraints 959 for r in list2_constraints: 960 r /= np.linalg.norm(r.ravel()) 961 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 968 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])) 976 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) 988 989 def __repr__(self): 990 constraints = repr(self.constraints) 991 return 'FixInternals(_copy_init=%s, epsilon=%s)' % (constraints, 992 repr(self.epsilon)) 993 994 def __str__(self): 995 return '\n'.join([repr(c) for c in self.constraints]) 996 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 1010 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 1020 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) 1026 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) 1030 1031 class FixBondCombo(FixInternalsBase): 1032 """Constraint subobject for fixing linear combination of bond lengths 1033 within FixInternals. 1034 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) 1042 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) 1051 1052 def __repr__(self): 1053 return 'FixBondCombo({}, {}, {})'.format(repr(self.targetvalue), 1054 self.indices, self.coefs) 1055 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) 1062 1063 def __repr__(self): 1064 return 'FixBondLengthAlt({}, {})'.format(self.targetvalue, 1065 *self.indices) 1066 1067 class FixAngleCombo(FixInternalsBase): 1068 """Constraint subobject for fixing linear combination of angles 1069 within FixInternals. 1070 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 1077 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) 1083 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) 1090 1091 def __repr__(self): 1092 return 'FixAngleCombo({}, {}, {})'.format(self.targetvalue, 1093 self.indices, self.coefs) 1094 1095 class FixAngle(FixAngleCombo): 1096 """Constraint object for fixing an angle within 1097 FixInternals using the SHAKE algorithm. 1098 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) 1106 1107 def __repr__(self): 1108 return 'FixAngle({}, {})'.format(self.targetvalue, *self.indices) 1109 1110 class FixDihedralCombo(FixInternalsBase): 1111 """Constraint subobject for fixing linear combination of dihedrals 1112 within FixInternals. 1113 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 1121 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) 1127 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) 1134 1135 def __repr__(self): 1136 return 'FixDihedralCombo({}, {}, {})'.format(self.targetvalue, 1137 self.indices, 1138 self.coefs) 1139 1140 class FixDihedral(FixDihedralCombo): 1141 """Constraint object for fixing a dihedral angle using 1142 the SHAKE algorithm. This one allows also other constraints. 1143 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) 1151 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) 1158 1159 def __repr__(self): 1160 return 'FixDihedral({}, {})'.format(self.targetvalue, *self.indices) 1161 1162 1163class FixParametricRelations(FixConstraint): 1164 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 1175 1176 These constraints are based off the work in: https://arxiv.org/abs/1908.01610 1177 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. 1182 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. 1186 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. 1191 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) 1208 1209 assert self.const_shift.shape[0] == 3*len(self.indices) 1210 assert self.Jacobian.shape[0] == 3*len(self.indices) 1211 1212 self.eps = eps 1213 self.use_cell = use_cell 1214 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] 1223 1224 self.params = params 1225 1226 self.Jacobian_inv = np.linalg.inv(self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T 1227 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 1231 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. 1242 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 ] 1251 1252 For diamond are: 1253 params = [] 1254 expressions = [ 1255 "0.0", "0.0", "0.0", 1256 "0.25", "0.25", "0.25", 1257 ], 1258 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 ] 1271 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 1281 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)) 1294 1295 for expr_ind, expression in enumerate(expressions): 1296 expression = expression.strip() 1297 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:] 1304 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}" 1307 1308 param_dct = dict() 1309 param_map = dict() 1310 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 1316 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]) 1321 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.") 1328 1329 const_shift[expr_ind] = float(eval_expression(expression, param_dct)) 1330 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 1340 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 1347 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) 1359 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) 1370 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 1376 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 += "-" 1385 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) 1390 1391 exp += param_exp 1392 1393 expressions.append(exp) 1394 return np.array(expressions).reshape((-1, 3)) 1395 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 } 1409 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]) 1416 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 = "[]" 1423 1424 return '{:s}({:s}, {:s}, ..., {:e})'.format( 1425 type(self).__name__, 1426 indices_str, 1427 params_str, 1428 self.eps 1429 ) 1430 1431 1432class FixScaledParametricRelations(FixParametricRelations): 1433 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 1443 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 ) 1454 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)) 1460 1461 return cell.cartesian_positions(scaled) 1462 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]) 1471 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) 1478 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 ) 1486 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 1493 1494 jacobian = cart2frac_jacob @ self.Jacobian 1495 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T 1496 1497 reduced_forces = jacobian.T @ forces.flatten() 1498 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3) 1499 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 1505 1506 1507class FixCartesianParametricRelations(FixParametricRelations): 1508 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 ) 1527 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 1533 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 ) 1542 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 ) 1551 1552 def adjust_forces(self, atoms, forces): 1553 """Adjust forces of the atoms to match the constraints""" 1554 if self.use_cell: 1555 return 1556 1557 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten() 1558 forces[self.indices] = (self.Jacobian_inv.T @ forces_reduced).reshape(-1, 3) 1559 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 ) 1568 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 1573 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) 1577 1578 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3) 1579 1580 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.""" 1584 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. 1591 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. 1605 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 """ 1614 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 1630 1631 def get_removed_dof(self, atoms): 1632 return 0 1633 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 1650 1651 def adjust_positions(self, atoms, newpositions): 1652 pass 1653 1654 def adjust_momenta(self, atoms, momenta): 1655 pass 1656 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 1685 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. 1711 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 1719 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 1740 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 1748 1749 1750class ExternalForce(FixConstraint): 1751 """Constraint object for pulling two atoms apart by an external force. 1752 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: 1756 1757 >>> con1 = ExternalForce(atom1, atom2, f_ext) 1758 >>> con2 = FixBondLength(atom3, atom4) 1759 >>> atoms.set_constraint([con1, con2]) 1760 1761 see ase/test/external_force.py""" 1762 1763 def __init__(self, a1, a2, f_ext): 1764 self.indices = [a1, a2] 1765 self.external_force = f_ext 1766 1767 def get_removed_dof(self, atoms): 1768 return 0 1769 1770 def adjust_positions(self, atoms, new): 1771 pass 1772 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) 1777 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 1781 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 1792 1793 def __repr__(self): 1794 return 'ExternalForce(%d, %d, %f)' % (self.indices[0], 1795 self.indices[1], 1796 self.external_force) 1797 1798 def todict(self): 1799 return {'name': 'ExternalForce', 1800 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1], 1801 'f_ext': self.external_force}} 1802 1803 1804class MirrorForce(FixConstraint): 1805 """Constraint object for mirroring the force between two atoms. 1806 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*. 1816 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. 1831 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: 1837 1838 >>> con1 = MirrorForce(atom1, atom2) 1839 >>> con2 = FixBondLength(atom3, atom4) 1840 >>> atoms.set_constraint([con1, con2]) 1841 1842 """ 1843 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 1849 1850 def adjust_positions(self, atoms, new): 1851 pass 1852 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) 1881 1882 def index_shuffle(self, atoms, ind): 1883 """Shuffle the indices of the two atoms in this constraint 1884 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 1894 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) 1899 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}} 1905 1906 1907class MirrorTorque(FixConstraint): 1908 """Constraint object for mirroring the torque acting on a dihedral 1909 angle defined by four atoms. 1910 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*. 1920 1921 This constraint can be used to find 1922 transition states of cis-trans-isomerization. 1923 1924 a1 a4 1925 | | 1926 a2 __ a3 1927 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. 1946 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: 1952 1953 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4) 1954 >>> con2 = FixBondLength(atom5, atom6) 1955 >>> atoms.set_constraint([con1, con2]) 1956 1957 """ 1958 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 1965 1966 def adjust_positions(self, atoms, new): 1967 pass 1968 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] 1979 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) 2018 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) 2028 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) 2033 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}} 2040 2041 2042class Filter: 2043 """Subset filter class.""" 2044 2045 def __init__(self, atoms, indices=None, mask=None): 2046 """Filter atoms. 2047 2048 This filter can be used to hide degrees of freedom in an Atoms 2049 object. 2050 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. 2058 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 """ 2063 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 2068 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".') 2073 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) 2080 2081 def iterimages(self): 2082 # Present the real atoms object to Trajectory and friends 2083 return self.atoms.iterimages() 2084 2085 def get_cell(self): 2086 """Returns the computational cell. 2087 2088 The computational cell is the same as for the original system. 2089 """ 2090 return self.atoms.get_cell() 2091 2092 def get_pbc(self): 2093 """Returns the periodic boundary conditions. 2094 2095 The boundary conditions are the same as for the original system. 2096 """ 2097 return self.atoms.get_pbc() 2098 2099 def get_positions(self): 2100 'Return the positions of the visible atoms.' 2101 return self.atoms.get_positions()[self.index] 2102 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) 2108 2109 positions = property(get_positions, set_positions, 2110 doc='Positions of the atoms') 2111 2112 def get_momenta(self): 2113 'Return the momenta of the visible atoms.' 2114 return self.atoms.get_momenta()[self.index] 2115 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) 2121 2122 def get_atomic_numbers(self): 2123 'Return the atomic numbers of the visible atoms.' 2124 return self.atoms.get_atomic_numbers()[self.index] 2125 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) 2131 2132 def get_tags(self): 2133 'Return the tags of the visible atoms.' 2134 return self.atoms.get_tags()[self.index] 2135 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) 2141 2142 def get_forces(self, *args, **kwargs): 2143 return self.atoms.get_forces(*args, **kwargs)[self.index] 2144 2145 def get_stress(self, *args, **kwargs): 2146 return self.atoms.get_stress(*args, **kwargs) 2147 2148 def get_stresses(self, *args, **kwargs): 2149 return self.atoms.get_stresses(*args, **kwargs)[self.index] 2150 2151 def get_masses(self): 2152 return self.atoms.get_masses()[self.index] 2153 2154 def get_potential_energy(self, **kwargs): 2155 """Calculate potential energy. 2156 2157 Returns the potential energy of the full system. 2158 """ 2159 return self.atoms.get_potential_energy(**kwargs) 2160 2161 def get_chemical_symbols(self): 2162 return self.atoms.get_chemical_symbols() 2163 2164 def get_initial_magnetic_moments(self): 2165 return self.atoms.get_initial_magnetic_moments() 2166 2167 def get_calculator(self): 2168 """Returns the calculator. 2169 2170 WARNING: The calculator is unaware of this filter, and sees a 2171 different number of atoms. 2172 """ 2173 return self.atoms.calc 2174 2175 @property 2176 def calc(self): 2177 return self.atoms.calc 2178 2179 def get_celldisp(self): 2180 return self.atoms.get_celldisp() 2181 2182 def has(self, name): 2183 'Check for existence of array.' 2184 return self.atoms.has(name) 2185 2186 def __len__(self): 2187 'Return the number of movable atoms.' 2188 return self.n 2189 2190 def __getitem__(self, i): 2191 'Return an atom.' 2192 return self.atoms[self.index[i]] 2193 2194 2195class StrainFilter(Filter): 2196 """Modify the supercell while keeping the scaled positions fixed. 2197 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. 2201 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. 2206 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. 2210 2211 """ 2212 2213 def __init__(self, atoms, mask=None, include_ideal_gas=False): 2214 """Create a filter applying a homogeneous strain to a list of atoms. 2215 2216 The first argument, atoms, is the atoms object. 2217 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]. 2222 2223 """ 2224 2225 self.strain = np.zeros(6) 2226 self.include_ideal_gas = include_ideal_gas 2227 2228 if mask is None: 2229 mask = np.ones(6) 2230 else: 2231 mask = np.array(mask) 2232 2233 Filter.__init__(self, atoms, mask=mask) 2234 self.mask = mask 2235 self.origcell = atoms.get_cell() 2236 2237 def get_positions(self): 2238 return self.strain.reshape((2, 3)).copy() 2239 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]]]) 2245 2246 self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True) 2247 self.strain[:] = new 2248 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)) 2252 2253 def has(self, x): 2254 return self.atoms.has(x) 2255 2256 def __len__(self): 2257 return 2 2258 2259 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. 2269 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. 2273 2274 - True = relax to zero 2275 - False = fixed, ignore this component 2276 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. 2281 2282 For full details see: 2283 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras, 2284 Phys. Rev. B 59, 235 (1999) 2285 2286 You can still use constraints on the atoms, e.g. FixAtoms, to control 2287 the relaxation of the atoms. 2288 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) 2293 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: 2297 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) 2304 2305 Helpful conversion table: 2306 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 2312 2313 Additional optional arguments: 2314 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. 2322 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). 2326 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). 2334 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 """ 2339 2340 Filter.__init__(self, atoms, indices=range(len(atoms))) 2341 self.atoms = atoms 2342 self.orig_cell = atoms.get_cell() 2343 self.stress = None 2344 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,)') 2354 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 2363 2364 def deform_grad(self): 2365 return np.linalg.solve(self.orig_cell, self.atoms.cell).T 2366 2367 def get_positions(self): 2368 """ 2369 this returns an array with shape (natoms + 3,3). 2370 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 """ 2375 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 2387 2388 def set_positions(self, new, **kwargs): 2389 """ 2390 new is an array with shape (natoms+3,3). 2391 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. 2394 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 """ 2399 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) 2415 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() 2423 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. 2428 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 """ 2433 2434 stress = self.atoms.get_stress(**kwargs) 2435 atoms_forces = self.atoms.get_forces(**kwargs) 2436 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 2443 2444 if self.hydrostatic_strain: 2445 vtr = virial.trace() 2446 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0]) 2447 2448 # Zero out components corresponding to fixed lattice elements 2449 if (self.mask != 1.0).any(): 2450 virial *= self.mask 2451 2452 if self.constant_volume: 2453 vtr = virial.trace() 2454 np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0) 2455 2456 natoms = len(self.atoms) 2457 forces = np.zeros((natoms + 3, 3)) 2458 forces[:natoms] = atoms_forces 2459 forces[natoms:] = virial / self.cell_factor 2460 2461 self.stress = -full_3x3_to_voigt_6_stress(virial) / volume 2462 return forces 2463 2464 def get_stress(self): 2465 raise PropertyNotImplementedError 2466 2467 def has(self, x): 2468 return self.atoms.has(x) 2469 2470 def __len__(self): 2471 return (len(self.atoms) + 3) 2472 2473 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. 2483 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. 2487 2488 - True = relax to zero 2489 - False = fixed, ignore this component 2490 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. 2495 2496 For full details see: 2497 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras, 2498 Phys. Rev. B 59, 235 (1999) 2499 2500 You can still use constraints on the atoms, e.g. FixAtoms, to control 2501 the relaxation of the atoms. 2502 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) 2507 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: 2511 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) 2518 2519 Helpful conversion table: 2520 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 2526 2527 Additional optional arguments: 2528 2529 cell_factor: (DEPRECATED) 2530 Retained for backwards compatibility, but no longer used. 2531 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). 2535 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. 2540 2541 scalar_pressure: float (default 0.0) 2542 Applied pressure to use for enthalpy pV term. As above, this 2543 breaks energy/force consistency. 2544 2545 Implementation details: 2546 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 2549 2550 We decompose the deformation gradient as 2551 2552 F = exp(U) F0 2553 x = F * F0^{-1} z = exp(U) z 2554 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). 2558 2559 \phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) . 2560 ( L(U, V) exp(-U) exp(U) z ) 2561 2562 where 2563 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') 2568 2569 where the : operator represents double contraction, 2570 i.e. A:B = trace(A'B), and 2571 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) 2583 2584 This means we can write 2585 2586 d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V 2587 2588 and therefore the contribution to the gradient of the energy is 2589 2590 \nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij 2591 """ 2592 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 2600 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 2606 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) 2612 2613 def get_forces(self, **kwargs): 2614 forces = UnitCellFilter.get_forces(self, **kwargs) 2615 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)) 2622 2623 cur_deform_grad = self.deform_grad() 2624 cur_deform_grad_log = logm(cur_deform_grad) 2625 2626 if self.hydrostatic_strain: 2627 vtr = virial.trace() 2628 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0]) 2629 2630 # Zero out components corresponding to fixed lattice elements 2631 if (self.mask != 1.0).any(): 2632 virial *= self.mask 2633 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 2645 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 2655 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) 2667 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 2673