Coverage for ase/constraints.py: 87.36%

1# fmt: off

3"""Constraints"""

4from __future__ import annotations

6from copy import deepcopy

7from typing import Any, Sequence

8from warnings import warn

10import numpy as np

12from ase import Atoms

13from ase.geometry import (

14 conditional_find_mic,

15 find_mic,

16 get_angles,

17 get_angles_derivatives,

18 get_dihedrals,

19 get_dihedrals_derivatives,

20 get_distances_derivatives,

21 wrap_positions,

22)

23from ase.spacegroup.symmetrize import (

24 prep_symmetry,

25 refine_symmetry,

26 symmetrize_rank1,

27 symmetrize_rank2,

28)

29from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress

30from ase.utils.parsemath import eval_expression

32__all__ = [

33 'FixCartesian', 'FixBondLength', 'FixedMode',

34 'FixAtoms', 'FixScaled', 'FixCom', 'FixSubsetCom', 'FixedPlane',

35 'FixConstraint', 'FixedLine', 'FixBondLengths', 'FixLinearTriatomic',

36 'FixInternals', 'Hookean', 'ExternalForce', 'MirrorForce', 'MirrorTorque',

37 'FixScaledParametricRelations', 'FixCartesianParametricRelations',

38 'FixSymmetry']

41def dict2constraint(dct: dict[str, Any]) -> FixConstraint:

42 """Convert dictionary to ASE `FixConstraint` object."""

43 if dct['name'] not in __all__:

44 raise ValueError

45 # address backward-compatibility breaking between ASE 3.22.0 and 3.23.0

46 # https://gitlab.com/ase/ase/-/merge_requests/3786

47 if dct['name'] in {'FixedLine', 'FixedPlane'} and 'a' in dct['kwargs']:

48 dct = deepcopy(dct)

49 dct['kwargs']['indices'] = dct['kwargs'].pop('a')

50 return globals()[dct['name']](**dct['kwargs'])

53def slice2enlist(s, n):

54 """Convert a slice object into a list of (new, old) tuples."""

55 if isinstance(s, slice):

56 return enumerate(range(*s.indices(n)))

57 return enumerate(s)

60def constrained_indices(atoms, only_include=None):

61 """Returns a list of indices for the atoms that are constrained

62 by a constraint that is applied. By setting only_include to a

63 specific type of constraint you can make it only look for that

64 given constraint.

65 """

66 indices = []

67 for constraint in atoms.constraints:

68 if only_include is not None:

69 if not isinstance(constraint, only_include):

70 continue

71 indices.extend(np.array(constraint.get_indices()))

72 return np.array(np.unique(indices))

75class FixConstraint:

76 """Base class for classes that fix one or more atoms in some way."""

78 def index_shuffle(self, atoms: Atoms, ind):

79 """Change the indices.

81 When the ordering of the atoms in the Atoms object changes,

82 this method can be called to shuffle the indices of the

83 constraints.

85 ind -- List or tuple of indices.

87 """

88 raise NotImplementedError

90 def repeat(self, m: int, n: int):

91 """ basic method to multiply by m, needs to know the length

92 of the underlying atoms object for the assignment of

93 multiplied constraints to work.

94 """

95 msg = ("Repeat is not compatible with your atoms' constraints."

96 ' Use atoms.set_constraint() before calling repeat to '

97 'remove your constraints.')

98 raise NotImplementedError(msg)

100 def get_removed_dof(self, atoms: Atoms):

101 """Get number of removed degrees of freedom due to constraint."""

102 raise NotImplementedError

103

104 def adjust_positions(self, atoms: Atoms, new):

105 """Adjust positions."""

106

107 def adjust_momenta(self, atoms: Atoms, momenta):

108 """Adjust momenta."""

109 # The default is in identical manner to forces.

110 # TODO: The default is however not always reasonable.

111 self.adjust_forces(atoms, momenta)

112

113 def adjust_forces(self, atoms: Atoms, forces):

114 """Adjust forces."""

115

116 def copy(self):

117 """Copy constraint."""

118 return dict2constraint(self.todict().copy())

119

120 def todict(self):

121 """Convert constraint to dictionary."""

122

123

124class IndexedConstraint(FixConstraint):

125 def __init__(self, indices=None, mask=None):

126 """Constrain chosen atoms.

127

128 Parameters

129 ----------

130 indices : sequence of int

131 Indices for those atoms that should be constrained.

132 mask : sequence of bool

133 One boolean per atom indicating if the atom should be

134 constrained or not.

135 """

136

137 if mask is not None:

138 if indices is not None:

139 raise ValueError('Use only one of "indices" and "mask".')

140 indices = mask

141 indices = np.atleast_1d(indices)

142 if np.ndim(indices) > 1:

143 raise ValueError('indices has wrong amount of dimensions. '

144 f'Got {np.ndim(indices)}, expected ndim <= 1')

145

146 if indices.dtype == bool:

147 indices = np.arange(len(indices))[indices]

148 elif len(indices) == 0:

149 indices = np.empty(0, dtype=int)

150 elif not np.issubdtype(indices.dtype, np.integer):

151 raise ValueError('Indices must be integers or boolean mask, '

152 f'not dtype={indices.dtype}')

153

154 if len(set(indices)) < len(indices):

155 raise ValueError(

156 'The indices array contains duplicates. '

157 'Perhaps you want to specify a mask instead, but '

158 'forgot the mask= keyword.')

159

160 self.index = indices

161

162 def index_shuffle(self, atoms, ind):

163 # See docstring of superclass

164 index = []

165

166 # Resolve negative indices:

167 actual_indices = set(np.arange(len(atoms))[self.index])

168

169 for new, old in slice2enlist(ind, len(atoms)):

170 if old in actual_indices:

171 index.append(new)

172 if len(index) == 0:

173 raise IndexError('All indices in FixAtoms not part of slice')

174 self.index = np.asarray(index, int)

175 # XXX make immutable

176

177 def get_indices(self):

178 return self.index.copy()

179

180 def repeat(self, m, n):

181 i0 = 0

182 natoms = 0

183 if isinstance(m, int):

184 m = (m, m, m)

185 index_new = []

186 for _ in range(m[2]):

187 for _ in range(m[1]):

188 for _ in range(m[0]):

189 i1 = i0 + n

190 index_new += [i + natoms for i in self.index]

191 i0 = i1

192 natoms += n

193 self.index = np.asarray(index_new, int)

194 # XXX make immutable

195 return self

196

197 def delete_atoms(self, indices, natoms):

198 """Removes atoms from the index array, if present.

199

200 Required for removing atoms with existing constraint.

201 """

202

203 i = np.zeros(natoms, int) - 1

204 new = np.delete(np.arange(natoms), indices)

205 i[new] = np.arange(len(new))

206 index = i[self.index]

207 self.index = index[index >= 0]

208 # XXX make immutable

209 if len(self.index) == 0:

210 return None

211 return self

212

213

214class FixAtoms(IndexedConstraint):

215 """Fix chosen atoms.

216

217 Examples

218 --------

219 Fix all Copper atoms:

220

221 >>> from ase.build import bulk

222

223 >>> atoms = bulk('Cu', 'fcc', a=3.6)

224 >>> mask = (atoms.symbols == 'Cu')

225 >>> c = FixAtoms(mask=mask)

226 >>> atoms.set_constraint(c)

227

228 Fix all atoms with z-coordinate less than 1.0 Angstrom:

229

230 >>> c = FixAtoms(mask=atoms.positions[:, 2] < 1.0)

231 >>> atoms.set_constraint(c)

232 """

233

234 def get_removed_dof(self, atoms):

235 return 3 * len(self.index)

236

237 def adjust_positions(self, atoms, new):

238 new[self.index] = atoms.positions[self.index]

239

240 def adjust_forces(self, atoms, forces):

241 forces[self.index] = 0.0

242

243 def __repr__(self):

244 clsname = type(self).__name__

245 indices = ints2string(self.index)

246 return f'{clsname}(indices={indices})'

247

248 def todict(self):

249 return {'name': 'FixAtoms',

250 'kwargs': {'indices': self.index.tolist()}}

251

252

253class FixCom(FixConstraint):

254 """Constraint class for fixing the center of mass."""

255

256 index = slice(None) # all atoms

257

258 def get_removed_dof(self, atoms):

259 return 3

260

261 def adjust_positions(self, atoms, new):

262 masses = atoms.get_masses()[self.index]

263 old_cm = atoms.get_center_of_mass(indices=self.index)

264 new_cm = masses @ new[self.index] / masses.sum()

265 diff = old_cm - new_cm

266 new += diff

267

268 def adjust_momenta(self, atoms, momenta):

269 """Adjust momenta so that the center-of-mass velocity is zero."""

270 masses = atoms.get_masses()[self.index]

271 velocity_com = momenta[self.index].sum(axis=0) / masses.sum()

272 momenta[self.index] -= masses[:, None] * velocity_com

273

274 def adjust_forces(self, atoms, forces):

275 # Eqs. (3) and (7) in https://doi.org/10.1021/jp9722824

276 masses = atoms.get_masses()[self.index]

277 lmd = masses @ forces[self.index] / sum(masses**2)

278 forces[self.index] -= masses[:, None] * lmd

279

280 def todict(self):

281 return {'name': 'FixCom',

282 'kwargs': {}}

283

284

285class FixSubsetCom(FixCom, IndexedConstraint):

286 """Constraint class for fixing the center of mass of a subset of atoms."""

287

288 def __init__(self, indices):

289 super().__init__(indices=indices)

290

291 def todict(self):

292 return {'name': self.__class__.__name__,

293 'kwargs': {'indices': self.index.tolist()}}

294

295

296def ints2string(x, threshold=None):

297 """Convert ndarray of ints to string."""

298 if threshold is None or len(x) <= threshold:

299 return str(x.tolist())

300 return str(x[:threshold].tolist())[:-1] + ', ...]'

301

302

303class FixBondLengths(FixConstraint):

304 maxiter = 500

305

306 def __init__(self, pairs, tolerance=1e-13,

307 bondlengths=None, iterations=None):

308 """iterations:

309 Ignored"""

310 self.pairs = np.asarray(pairs)

311 self.tolerance = tolerance

312 self.bondlengths = bondlengths

313

314 def get_removed_dof(self, atoms):

315 return len(self.pairs)

316

317 def adjust_positions(self, atoms, new):

318 old = atoms.positions

319 masses = atoms.get_masses()

320

321 if self.bondlengths is None:

322 self.bondlengths = self.initialize_bond_lengths(atoms)

323

324 for i in range(self.maxiter):

325 converged = True

326 for j, ab in enumerate(self.pairs):

327 a = ab[0]

328 b = ab[1]

329 cd = self.bondlengths[j]

330 r0 = old[a] - old[b]

331 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)

332 d1 = new[a] - new[b] - r0 + d0

333 m = 1 / (1 / masses[a] + 1 / masses[b])

334 x = 0.5 * (cd**2 - np.dot(d1, d1)) / np.dot(d0, d1)

335 if abs(x) > self.tolerance:

336 new[a] += x * m / masses[a] * d0

337 new[b] -= x * m / masses[b] * d0

338 converged = False

339 if converged:

340 break

341 else:

342 raise RuntimeError('Did not converge')

343

344 def adjust_momenta(self, atoms, p):

345 old = atoms.positions

346 masses = atoms.get_masses()

347

348 if self.bondlengths is None:

349 self.bondlengths = self.initialize_bond_lengths(atoms)

350

351 for i in range(self.maxiter):

352 converged = True

353 for j, ab in enumerate(self.pairs):

354 a = ab[0]

355 b = ab[1]

356 cd = self.bondlengths[j]

357 d = old[a] - old[b]

358 d, _ = find_mic(d, atoms.cell, atoms.pbc)

359 dv = p[a] / masses[a] - p[b] / masses[b]

360 m = 1 / (1 / masses[a] + 1 / masses[b])

361 x = -np.dot(dv, d) / cd**2

362 if abs(x) > self.tolerance:

363 p[a] += x * m * d

364 p[b] -= x * m * d

365 converged = False

366 if converged:

367 break

368 else:

369 raise RuntimeError('Did not converge')

370

371 def adjust_forces(self, atoms, forces):

372 self.constraint_forces = -forces

373 self.adjust_momenta(atoms, forces)

374 self.constraint_forces += forces

375

376 def initialize_bond_lengths(self, atoms):

377 bondlengths = np.zeros(len(self.pairs))

378

379 for i, ab in enumerate(self.pairs):

380 bondlengths[i] = atoms.get_distance(ab[0], ab[1], mic=True)

381

382 return bondlengths

383

384 def get_indices(self):

385 return np.unique(self.pairs.ravel())

386

387 def todict(self):

388 return {'name': 'FixBondLengths',

389 'kwargs': {'pairs': self.pairs.tolist(),

390 'tolerance': self.tolerance}}

391

392 def index_shuffle(self, atoms, ind):

393 """Shuffle the indices of the two atoms in this constraint"""

394 map = np.zeros(len(atoms), int)

395 map[ind] = 1

396 n = map.sum()

397 map[:] = -1

398 map[ind] = range(n)

399 pairs = map[self.pairs]

400 self.pairs = pairs[(pairs != -1).all(1)]

401 if len(self.pairs) == 0:

402 raise IndexError('Constraint not part of slice')

403

404

405def FixBondLength(a1, a2):

406 """Fix distance between atoms with indices a1 and a2."""

407 return FixBondLengths([(a1, a2)])

408

409

410class FixLinearTriatomic(FixConstraint):

411 """Holonomic constraints for rigid linear triatomic molecules."""

412

413 def __init__(self, triples):

414 """Apply RATTLE-type bond constraints between outer atoms n and m

415 and linear vectorial constraints to the position of central

416 atoms o to fix the geometry of linear triatomic molecules of the

417 type:

418

419 n--o--m

420

421 Parameters:

422

423 triples: list

424 Indices of the atoms forming the linear molecules to constrain

425 as triples. Sequence should be (n, o, m) or (m, o, n).

426

427 When using these constraints in molecular dynamics or structure

428 optimizations, atomic forces need to be redistributed within a

429 triple. The function redistribute_forces_optimization implements

430 the redistribution of forces for structure optimization, while

431 the function redistribute_forces_md implements the redistribution

432 for molecular dynamics.

433

434 References:

435

436 Ciccotti et al. Molecular Physics 47 (1982)

437 :doi:`10.1080/00268978200100942`

438 """

439 self.triples = np.asarray(triples)

440 if self.triples.shape[1] != 3:

441 raise ValueError('"triples" has wrong size')

442 self.bondlengths = None

443

444 def get_removed_dof(self, atoms):

445 return 4 * len(self.triples)

446

447 @property

448 def n_ind(self):

449 return self.triples[:, 0]

450

451 @property

452 def m_ind(self):

453 return self.triples[:, 2]

454

455 @property

456 def o_ind(self):

457 return self.triples[:, 1]

458

459 def initialize(self, atoms):

460 masses = atoms.get_masses()

461 self.mass_n, self.mass_m, self.mass_o = self.get_slices(masses)

462

463 self.bondlengths = self.initialize_bond_lengths(atoms)

464 self.bondlengths_nm = self.bondlengths.sum(axis=1)

465

466 C1 = self.bondlengths[:, ::-1] / self.bondlengths_nm[:, None]

467 C2 = (C1[:, 0] ** 2 * self.mass_o * self.mass_m +

468 C1[:, 1] ** 2 * self.mass_n * self.mass_o +

469 self.mass_n * self.mass_m)

470 C2 = C1 / C2[:, None]

471 C3 = self.mass_n * C1[:, 1] - self.mass_m * C1[:, 0]

472 C3 = C2 * self.mass_o[:, None] * C3[:, None]

473 C3[:, 1] *= -1

474 C3 = (C3 + 1) / np.vstack((self.mass_n, self.mass_m)).T

475 C4 = (C1[:, 0]**2 + C1[:, 1]**2 + 1)

476 C4 = C1 / C4[:, None]

477

478 self.C1 = C1

479 self.C2 = C2

480 self.C3 = C3

481 self.C4 = C4

482

483 def adjust_positions(self, atoms, new):

484 old = atoms.positions

485 new_n, new_m, new_o = self.get_slices(new)

486

487 if self.bondlengths is None:

488 self.initialize(atoms)

489

490 r0 = old[self.n_ind] - old[self.m_ind]

491 d0, _ = find_mic(r0, atoms.cell, atoms.pbc)

492 d1 = new_n - new_m - r0 + d0

493 a = np.einsum('ij,ij->i', d0, d0)

494 b = np.einsum('ij,ij->i', d1, d0)

495 c = np.einsum('ij,ij->i', d1, d1) - self.bondlengths_nm ** 2

496 g = (b - (b**2 - a * c)**0.5) / (a * self.C3.sum(axis=1))

497 g = g[:, None] * self.C3

498 new_n -= g[:, 0, None] * d0

499 new_m += g[:, 1, None] * d0

500 if np.allclose(d0, r0):

501 new_o = (self.C1[:, 0, None] * new_n

502 + self.C1[:, 1, None] * new_m)

503 else:

504 v1, _ = find_mic(new_n, atoms.cell, atoms.pbc)

505 v2, _ = find_mic(new_m, atoms.cell, atoms.pbc)

506 rb = self.C1[:, 0, None] * v1 + self.C1[:, 1, None] * v2

507 new_o = wrap_positions(rb, atoms.cell, atoms.pbc)

508

509 self.set_slices(new_n, new_m, new_o, new)

510

511 def adjust_momenta(self, atoms, p):

512 old = atoms.positions

513 p_n, p_m, p_o = self.get_slices(p)

514

515 if self.bondlengths is None:

516 self.initialize(atoms)

517

518 mass_nn = self.mass_n[:, None]

519 mass_mm = self.mass_m[:, None]

520 mass_oo = self.mass_o[:, None]

521

522 d = old[self.n_ind] - old[self.m_ind]

523 d, _ = find_mic(d, atoms.cell, atoms.pbc)

524 dv = p_n / mass_nn - p_m / mass_mm

525 k = np.einsum('ij,ij->i', dv, d) / self.bondlengths_nm ** 2

526 k = self.C3 / (self.C3.sum(axis=1)[:, None]) * k[:, None]

527 p_n -= k[:, 0, None] * mass_nn * d

528 p_m += k[:, 1, None] * mass_mm * d

529 p_o = (mass_oo * (self.C1[:, 0, None] * p_n / mass_nn +

530 self.C1[:, 1, None] * p_m / mass_mm))

531

532 self.set_slices(p_n, p_m, p_o, p)

533

534 def adjust_forces(self, atoms, forces):

535

536 if self.bondlengths is None:

537 self.initialize(atoms)

538

539 A = self.C4 * np.diff(self.C1)

540 A[:, 0] *= -1

541 A -= 1

542 B = np.diff(self.C4) / (A.sum(axis=1))[:, None]

543 A /= (A.sum(axis=1))[:, None]

544

545 self.constraint_forces = -forces

546 old = atoms.positions

547

548 fr_n, fr_m, fr_o = self.redistribute_forces_optimization(forces)

549

550 d = old[self.n_ind] - old[self.m_ind]

551 d, _ = find_mic(d, atoms.cell, atoms.pbc)

552 df = fr_n - fr_m

553 k = -np.einsum('ij,ij->i', df, d) / self.bondlengths_nm ** 2

554 forces[self.n_ind] = fr_n + k[:, None] * d * A[:, 0, None]

555 forces[self.m_ind] = fr_m - k[:, None] * d * A[:, 1, None]

556 forces[self.o_ind] = fr_o + k[:, None] * d * B

557

558 self.constraint_forces += forces

559

560 def redistribute_forces_optimization(self, forces):

561 """Redistribute forces within a triple when performing structure

562 optimizations.

563

564 The redistributed forces needs to be further adjusted using the

565 appropriate Lagrange multipliers as implemented in adjust_forces."""

566 forces_n, forces_m, forces_o = self.get_slices(forces)

567 C1_1 = self.C1[:, 0, None]

568 C1_2 = self.C1[:, 1, None]

569 C4_1 = self.C4[:, 0, None]

570 C4_2 = self.C4[:, 1, None]

571

572 fr_n = ((1 - C4_1 * C1_1) * forces_n -

573 C4_1 * (C1_2 * forces_m - forces_o))

574 fr_m = ((1 - C4_2 * C1_2) * forces_m -

575 C4_2 * (C1_1 * forces_n - forces_o))

576 fr_o = ((1 - 1 / (C1_1**2 + C1_2**2 + 1)) * forces_o +

577 C4_1 * forces_n + C4_2 * forces_m)

578

579 return fr_n, fr_m, fr_o

580

581 def redistribute_forces_md(self, atoms, forces, rand=False):

582 """Redistribute forces within a triple when performing molecular

583 dynamics.

584

585 When rand=True, use the equations for random force terms, as

586 used e.g. by Langevin dynamics, otherwise apply the standard

587 equations for deterministic forces (see Ciccotti et al. Molecular

588 Physics 47 (1982))."""

589 if self.bondlengths is None:

590 self.initialize(atoms)

591 forces_n, forces_m, forces_o = self.get_slices(forces)

592 C1_1 = self.C1[:, 0, None]

593 C1_2 = self.C1[:, 1, None]

594 C2_1 = self.C2[:, 0, None]

595 C2_2 = self.C2[:, 1, None]

596 mass_nn = self.mass_n[:, None]

597 mass_mm = self.mass_m[:, None]

598 mass_oo = self.mass_o[:, None]

599 if rand:

600 mr1 = (mass_mm / mass_nn) ** 0.5

601 mr2 = (mass_oo / mass_nn) ** 0.5

602 mr3 = (mass_nn / mass_mm) ** 0.5

603 mr4 = (mass_oo / mass_mm) ** 0.5

604 else:

605 mr1 = 1.0

606 mr2 = 1.0

607 mr3 = 1.0

608 mr4 = 1.0

609

610 fr_n = ((1 - C1_1 * C2_1 * mass_oo * mass_mm) * forces_n -

611 C2_1 * (C1_2 * mr1 * mass_oo * mass_nn * forces_m -

612 mr2 * mass_mm * mass_nn * forces_o))

613

614 fr_m = ((1 - C1_2 * C2_2 * mass_oo * mass_nn) * forces_m -

615 C2_2 * (C1_1 * mr3 * mass_oo * mass_mm * forces_n -

616 mr4 * mass_mm * mass_nn * forces_o))

617

618 self.set_slices(fr_n, fr_m, 0.0, forces)

619

620 def get_slices(self, a):

621 a_n = a[self.n_ind]

622 a_m = a[self.m_ind]

623 a_o = a[self.o_ind]

624

625 return a_n, a_m, a_o

626

627 def set_slices(self, a_n, a_m, a_o, a):

628 a[self.n_ind] = a_n

629 a[self.m_ind] = a_m

630 a[self.o_ind] = a_o

631

632 def initialize_bond_lengths(self, atoms):

633 bondlengths = np.zeros((len(self.triples), 2))

634

635 for i in range(len(self.triples)):

636 bondlengths[i, 0] = atoms.get_distance(self.n_ind[i],

637 self.o_ind[i], mic=True)

638 bondlengths[i, 1] = atoms.get_distance(self.o_ind[i],

639 self.m_ind[i], mic=True)

640

641 return bondlengths

642

643 def get_indices(self):

644 return np.unique(self.triples.ravel())

645

646 def todict(self):

647 return {'name': 'FixLinearTriatomic',

648 'kwargs': {'triples': self.triples.tolist()}}

649

650 def index_shuffle(self, atoms, ind):

651 """Shuffle the indices of the three atoms in this constraint"""

652 map = np.zeros(len(atoms), int)

653 map[ind] = 1

654 n = map.sum()

655 map[:] = -1

656 map[ind] = range(n)

657 triples = map[self.triples]

658 self.triples = triples[(triples != -1).all(1)]

659 if len(self.triples) == 0:

660 raise IndexError('Constraint not part of slice')

661

662

663class FixedMode(FixConstraint):

664 """Constrain atoms to move along directions orthogonal to

665 a given mode only. Initialize with a mode, such as one produced by

666 ase.vibrations.Vibrations.get_mode()."""

667

668 def __init__(self, mode):

669 mode = np.asarray(mode)

670 self.mode = (mode / np.sqrt((mode**2).sum())).reshape(-1)

671

672 def get_removed_dof(self, atoms):

673 return len(atoms)

674

675 def adjust_positions(self, atoms, newpositions):

676 newpositions = newpositions.ravel()

677 oldpositions = atoms.positions.ravel()

678 step = newpositions - oldpositions

679 newpositions -= self.mode * np.dot(step, self.mode)

680

681 def adjust_forces(self, atoms, forces):

682 forces = forces.ravel()

683 forces -= self.mode * np.dot(forces, self.mode)

684

685 def index_shuffle(self, atoms, ind):

686 eps = 1e-12

687 mode = self.mode.reshape(-1, 3)

688 excluded = np.ones(len(mode), dtype=bool)

689 excluded[ind] = False

690 if (abs(mode[excluded]) > eps).any():

691 raise IndexError('All nonzero parts of mode not in slice')

692 self.mode = mode[ind].ravel()

693

694 def get_indices(self):

695 # This function will never properly work because it works on all

696 # atoms and it has no idea how to tell how many atoms it is

697 # attached to. If it is being used, surely the user knows

698 # everything is being constrained.

699 return []

700

701 def todict(self):

702 return {'name': 'FixedMode',

703 'kwargs': {'mode': self.mode.tolist()}}

704

705 def __repr__(self):

706 return f'FixedMode({self.mode.tolist()})'

707

708

709def _normalize(direction):

710 if np.shape(direction) != (3,):

711 raise ValueError("len(direction) is {len(direction)}. Has to be 3")

712

713 direction = np.asarray(direction) / np.linalg.norm(direction)

714 return direction

715

716

717class FixedPlane(IndexedConstraint):

718 """

719 Constraint object for fixing chosen atoms to only move in a plane.

720

721 The plane is defined by its normal vector *direction*

722 """

723

724 def __init__(self, indices, direction):

725 """Constrain chosen atoms.

726

727 Parameters

728 ----------

729 indices : int or list of int

730 Index or indices for atoms that should be constrained

731 direction : list of 3 int

732 Direction of the normal vector

733

734 Examples

735 --------

736 Fix all Copper atoms to only move in the yz-plane:

737

738 >>> from ase.build import bulk

739 >>> from ase.constraints import FixedPlane

740

741 >>> atoms = bulk('Cu', 'fcc', a=3.6)

742 >>> c = FixedPlane(

743 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],

744 ... direction=[1, 0, 0],

745 ... )

746 >>> atoms.set_constraint(c)

747

748 or constrain a single atom with the index 0 to move in the xy-plane:

749

750 >>> c = FixedPlane(indices=0, direction=[0, 0, 1])

751 >>> atoms.set_constraint(c)

752 """

753 super().__init__(indices=indices)

754 self.dir = _normalize(direction)

755

756 def adjust_positions(self, atoms, newpositions):

757 step = newpositions[self.index] - atoms.positions[self.index]

758 newpositions[self.index] -= _projection(step, self.dir)

759

760 def adjust_forces(self, atoms, forces):

761 forces[self.index] -= _projection(forces[self.index], self.dir)

762

763 def get_removed_dof(self, atoms):

764 return len(self.index)

765

766 def todict(self):

767 return {

768 'name': 'FixedPlane',

769 'kwargs': {'indices': self.index.tolist(),

770 'direction': self.dir.tolist()}

771 }

772

773 def __repr__(self):

774 return f'FixedPlane(indices={self.index}, {self.dir.tolist()})'

775

776

777def _projection(vectors, direction):

778 dotprods = vectors @ direction

779 projection = direction[None, :] * dotprods[:, None]

780 return projection

781

782

783class FixedLine(IndexedConstraint):

784 """

785 Constrain an atom index or a list of atom indices to move on a line only.

786

787 The line is defined by its vector *direction*

788 """

789

790 def __init__(self, indices, direction):

791 """Constrain chosen atoms.

792

793 Parameters

794 ----------

795 indices : int or list of int

796 Index or indices for atoms that should be constrained

797 direction : list of 3 int

798 Direction of the vector defining the line

799

800 Examples

801 --------

802 Fix all Copper atoms to only move in the x-direction:

803

804 >>> from ase.constraints import FixedLine

805 >>> c = FixedLine(

806 ... indices=[atom.index for atom in atoms if atom.symbol == 'Cu'],

807 ... direction=[1, 0, 0],

808 ... )

809 >>> atoms.set_constraint(c)

810

811 or constrain a single atom with the index 0 to move in the z-direction:

812

813 >>> c = FixedLine(indices=0, direction=[0, 0, 1])

814 >>> atoms.set_constraint(c)

815 """

816 super().__init__(indices)

817 self.dir = _normalize(direction)

818

819 def adjust_positions(self, atoms, newpositions):

820 step = newpositions[self.index] - atoms.positions[self.index]

821 projection = _projection(step, self.dir)

822 newpositions[self.index] = atoms.positions[self.index] + projection

823

824 def adjust_forces(self, atoms, forces):

825 forces[self.index] = _projection(forces[self.index], self.dir)

826

827 def get_removed_dof(self, atoms):

828 return 2 * len(self.index)

829

830 def __repr__(self):

831 return f'FixedLine(indices={self.index}, {self.dir.tolist()})'

832

833 def todict(self):

834 return {

835 'name': 'FixedLine',

836 'kwargs': {'indices': self.index.tolist(),

837 'direction': self.dir.tolist()}

838 }

839

840

841class FixCartesian(IndexedConstraint):

842 """Fix atoms in the directions of the cartesian coordinates.

843

844 Parameters

845 ----------

846 a : Sequence[int]

847 Indices of atoms to be fixed.

848 mask : tuple[bool, bool, bool], default: (True, True, True)

849 Cartesian directions to be fixed. (False: unfixed, True: fixed)

850 """

851

852 def __init__(self, a, mask=(True, True, True)):

853 super().__init__(indices=a)

854 self.mask = np.asarray(mask, bool)

855

856 def get_removed_dof(self, atoms: Atoms):

857 return self.mask.sum() * len(self.index)

858

859 def adjust_positions(self, atoms: Atoms, new):

860 new[self.index] = np.where(

861 self.mask[None, :],

862 atoms.positions[self.index],

863 new[self.index],

864 )

865

866 def adjust_forces(self, atoms: Atoms, forces):

867 forces[self.index] *= ~self.mask[None, :]

868

869 def todict(self):

870 return {'name': 'FixCartesian',

871 'kwargs': {'a': self.index.tolist(),

872 'mask': self.mask.tolist()}}

873

874 def __repr__(self):

875 name = type(self).__name__

876 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'

877

878

879class FixScaled(IndexedConstraint):

880 """Fix atoms in the directions of the unit vectors.

881

882 Parameters

883 ----------

884 a : Sequence[int]

885 Indices of atoms to be fixed.

886 mask : tuple[bool, bool, bool], default: (True, True, True)

887 Cell directions to be fixed. (False: unfixed, True: fixed)

888 """

889

890 def __init__(self, a, mask=(True, True, True), cell=None):

891 # XXX The unused cell keyword is there for compatibility

892 # with old trajectory files.

893 super().__init__(indices=a)

894 self.mask = np.asarray(mask, bool)

895

896 def get_removed_dof(self, atoms: Atoms):

897 return self.mask.sum() * len(self.index)

898

899 def adjust_positions(self, atoms: Atoms, new):

900 cell = atoms.cell

901 scaled_old = cell.scaled_positions(atoms.positions[self.index])

902 scaled_new = cell.scaled_positions(new[self.index])

903 scaled_new[:, self.mask] = scaled_old[:, self.mask]

904 new[self.index] = cell.cartesian_positions(scaled_new)

905

906 def adjust_forces(self, atoms: Atoms, forces):

907 # Forces are covariant to the coordinate transformation,

908 # use the inverse transformations

909 cell = atoms.cell

910 scaled_forces = cell.cartesian_positions(forces[self.index])

911 scaled_forces *= -(self.mask - 1)

912 forces[self.index] = cell.scaled_positions(scaled_forces)

913

914 def todict(self):

915 return {'name': 'FixScaled',

916 'kwargs': {'a': self.index.tolist(),

917 'mask': self.mask.tolist()}}

918

919 def __repr__(self):

920 name = type(self).__name__

921 return f'{name}(indices={self.index.tolist()}, {self.mask.tolist()})'

922

923

924# TODO: Better interface might be to use dictionaries in place of very

925# nested lists/tuples

926class FixInternals(FixConstraint):

927 """Constraint object for fixing multiple internal coordinates.

928

929 Allows fixing bonds, angles, dihedrals as well as linear combinations

930 of bonds (bondcombos).

931

932 Please provide angular units in degrees using `angles_deg` and

933 `dihedrals_deg`.

934 Fixing planar angles is not supported at the moment.

935 """

936

937 def __init__(self, bonds=None, angles=None, dihedrals=None,

938 angles_deg=None, dihedrals_deg=None,

939 bondcombos=None,

940 mic=False, epsilon=1.e-7):

941 """

942 A constrained internal coordinate is defined as a nested list:

943 '[value, [atom indices]]'. The constraint is initialized with a list of

944 constrained internal coordinates, i.e. '[[value, [atom indices]], ...]'.

945 If 'value' is None, the current value of the coordinate is constrained.

946

947 Parameters

948 ----------

949 bonds: nested python list, optional

950 List with targetvalue and atom indices defining the fixed bonds,

951 i.e. [[targetvalue, [index0, index1]], ...]

952

953 angles_deg: nested python list, optional

954 List with targetvalue and atom indices defining the fixedangles,

955 i.e. [[targetvalue, [index0, index1, index3]], ...]

956

957 dihedrals_deg: nested python list, optional

958 List with targetvalue and atom indices defining the fixed dihedrals,

959 i.e. [[targetvalue, [index0, index1, index3]], ...]

960

961 bondcombos: nested python list, optional

962 List with targetvalue, atom indices and linear coefficient defining

963 the fixed linear combination of bonds,

964 i.e. [[targetvalue, [[index0, index1, coefficient_for_bond],

965 [index1, index2, coefficient_for_bond]]], ...]

966

967 mic: bool, optional, default: False

968 Minimum image convention.

969

970 epsilon: float, optional, default: 1e-7

971 Convergence criterion.

972 """

973 warn_msg = 'Please specify {} in degrees using the {} argument.'

974 if angles:

975 warn(warn_msg.format('angles', 'angle_deg'), FutureWarning)

976 angles = np.asarray(angles)

977 angles[:, 0] = angles[:, 0] / np.pi * 180

978 angles = angles.tolist()

979 else:

980 angles = angles_deg

981 if dihedrals:

982 warn(warn_msg.format('dihedrals', 'dihedrals_deg'), FutureWarning)

983 dihedrals = np.asarray(dihedrals)

984 dihedrals[:, 0] = dihedrals[:, 0] / np.pi * 180

985 dihedrals = dihedrals.tolist()

986 else:

987 dihedrals = dihedrals_deg

988

989 self.bonds = bonds or []

990 self.angles = angles or []

991 self.dihedrals = dihedrals or []

992 self.bondcombos = bondcombos or []

993 self.mic = mic

994 self.epsilon = epsilon

995

996 self.n = (len(self.bonds) + len(self.angles) + len(self.dihedrals)

997 + len(self.bondcombos))

998

999 # Initialize these at run-time:

1000 self.constraints = []

1001 self.initialized = False

1002

1003 def get_removed_dof(self, atoms):

1004 return self.n

1005

1006 def initialize(self, atoms):

1007 if self.initialized:

1008 return

1009 masses = np.repeat(atoms.get_masses(), 3)

1010 cell = None

1011 pbc = None

1012 if self.mic:

1013 cell = atoms.cell

1014 pbc = atoms.pbc

1015 self.constraints = []

1016 for data, ConstrClass in [(self.bonds, self.FixBondLengthAlt),

1017 (self.angles, self.FixAngle),

1018 (self.dihedrals, self.FixDihedral),

1019 (self.bondcombos, self.FixBondCombo)]:

1020 for datum in data:

1021 targetvalue = datum[0]

1022 if targetvalue is None: # set to current value

1023 targetvalue = ConstrClass.get_value(atoms, datum[1],

1024 self.mic)

1025 constr = ConstrClass(targetvalue, datum[1], masses, cell, pbc)

1026 self.constraints.append(constr)

1027 self.initialized = True

1028

1029 @staticmethod

1030 def get_bondcombo(atoms, indices, mic=False):

1031 """Convenience function to return the value of the bondcombo coordinate

1032 (linear combination of bond lengths) for the given Atoms object 'atoms'.

1033 Example: Get the current value of the linear combination of two bond

1034 lengths defined as `bondcombo = [[0, 1, 1.0], [2, 3, -1.0]]`."""

1035 c = sum(df[2] * atoms.get_distance(*df[:2], mic=mic) for df in indices)

1036 return c

1037

1038 def get_subconstraint(self, atoms, definition):

1039 """Get pointer to a specific subconstraint.

1040 Identification by its definition via indices (and coefficients)."""

1041 self.initialize(atoms)

1042 for subconstr in self.constraints:

1043 if isinstance(definition[0], Sequence): # Combo constraint

1044 defin = [d + [c] for d, c in zip(subconstr.indices,

1045 subconstr.coefs)]

1046 if defin == definition:

1047 return subconstr

1048 else: # identify primitive constraints by their indices

1049 if subconstr.indices == [definition]:

1050 return subconstr

1051 raise ValueError('Given `definition` not found on Atoms object.')

1052

1053 def shuffle_definitions(self, shuffle_dic, internal_type):

1054 dfns = [] # definitions

1055 for dfn in internal_type: # e.g. for bond in self.bonds

1056 append = True

1057 new_dfn = [dfn[0], list(dfn[1])]

1058 for old in dfn[1]:

1059 if old in shuffle_dic:

1060 new_dfn[1][dfn[1].index(old)] = shuffle_dic[old]

1061 else:

1062 append = False

1063 break

1064 if append:

1065 dfns.append(new_dfn)

1066 return dfns

1067

1068 def shuffle_combos(self, shuffle_dic, internal_type):

1069 dfns = [] # definitions

1070 for dfn in internal_type: # i.e. for bondcombo in self.bondcombos

1071 append = True

1072 all_indices = [idx[0:-1] for idx in dfn[1]]

1073 new_dfn = [dfn[0], list(dfn[1])]

1074 for i, indices in enumerate(all_indices):

1075 for old in indices:

1076 if old in shuffle_dic:

1077 new_dfn[1][i][indices.index(old)] = shuffle_dic[old]

1078 else:

1079 append = False

1080 break

1081 if not append:

1082 break

1083 if append:

1084 dfns.append(new_dfn)

1085 return dfns

1086

1087 def index_shuffle(self, atoms, ind):

1088 # See docstring of superclass

1089 self.initialize(atoms)

1090 shuffle_dic = dict(slice2enlist(ind, len(atoms)))

1091 shuffle_dic = {old: new for new, old in shuffle_dic.items()}

1092 self.bonds = self.shuffle_definitions(shuffle_dic, self.bonds)

1093 self.angles = self.shuffle_definitions(shuffle_dic, self.angles)

1094 self.dihedrals = self.shuffle_definitions(shuffle_dic, self.dihedrals)

1095 self.bondcombos = self.shuffle_combos(shuffle_dic, self.bondcombos)

1096 self.initialized = False

1097 self.initialize(atoms)

1098 if len(self.constraints) == 0:

1099 raise IndexError('Constraint not part of slice')

1100

1101 def get_indices(self):

1102 cons = []

1103 for dfn in self.bonds + self.dihedrals + self.angles:

1104 cons.extend(dfn[1])

1105 for dfn in self.bondcombos:

1106 for partial_dfn in dfn[1]:

1107 cons.extend(partial_dfn[0:-1]) # last index is the coefficient

1108 return list(set(cons))

1109

1110 def todict(self):

1111 return {'name': 'FixInternals',

1112 'kwargs': {'bonds': self.bonds,

1113 'angles_deg': self.angles,

1114 'dihedrals_deg': self.dihedrals,

1115 'bondcombos': self.bondcombos,

1116 'mic': self.mic,

1117 'epsilon': self.epsilon}}

1118

1119 def adjust_positions(self, atoms, newpos):

1120 self.initialize(atoms)

1121 for constraint in self.constraints:

1122 constraint.setup_jacobian(atoms.positions)

1123 for _ in range(50):

1124 maxerr = 0.0

1125 for constraint in self.constraints:

1126 constraint.adjust_positions(atoms.positions, newpos)

1127 maxerr = max(abs(constraint.sigma), maxerr)

1128 if maxerr < self.epsilon:

1129 return

1130 msg = 'FixInternals.adjust_positions did not converge.'

1131 if any(constr.targetvalue > 175. or constr.targetvalue < 5. for constr

1132 in self.constraints if isinstance(constr, self.FixAngle)):

1133 msg += (' This may be caused by an almost planar angle.'

1134 ' Support for planar angles would require the'

1135 ' implementation of ghost, i.e. dummy, atoms.'

1136 ' See issue #868.')

1137 raise ValueError(msg)

1138

1139 def adjust_forces(self, atoms, forces):

1140 """Project out translations and rotations and all other constraints"""

1141 self.initialize(atoms)

1142 positions = atoms.positions

1143 N = len(forces)

1144 list2_constraints = list(np.zeros((6, N, 3)))

1145 tx, ty, tz, rx, ry, rz = list2_constraints

1146

1147 list_constraints = [r.ravel() for r in list2_constraints]

1148

1149 tx[:, 0] = 1.0

1150 ty[:, 1] = 1.0

1151 tz[:, 2] = 1.0

1152 ff = forces.ravel()

1153

1154 # Calculate the center of mass

1155 center = positions.sum(axis=0) / N

1156

1157 rx[:, 1] = -(positions[:, 2] - center[2])

1158 rx[:, 2] = positions[:, 1] - center[1]

1159 ry[:, 0] = positions[:, 2] - center[2]

1160 ry[:, 2] = -(positions[:, 0] - center[0])

1161 rz[:, 0] = -(positions[:, 1] - center[1])

1162 rz[:, 1] = positions[:, 0] - center[0]

1163

1164 # Normalizing transl., rotat. constraints

1165 for r in list2_constraints:

1166 r /= np.linalg.norm(r.ravel())

1167

1168 # Add all angle, etc. constraint vectors

1169 for constraint in self.constraints:

1170 constraint.setup_jacobian(positions)

1171 constraint.adjust_forces(positions, forces)

1172 list_constraints.insert(0, constraint.jacobian)

1173 # QR DECOMPOSITION - GRAM SCHMIDT

1174

1175 list_constraints = [r.ravel() for r in list_constraints]

1176 aa = np.column_stack(list_constraints)

1177 (aa, _bb) = np.linalg.qr(aa)

1178 # Projection

1179 hh = []

1180 for i, constraint in enumerate(self.constraints):

1181 hh.append(aa[:, i] * np.vstack(aa[:, i]))

1182

1183 txx = aa[:, self.n] * np.vstack(aa[:, self.n])

1184 tyy = aa[:, self.n + 1] * np.vstack(aa[:, self.n + 1])

1185 tzz = aa[:, self.n + 2] * np.vstack(aa[:, self.n + 2])

1186 rxx = aa[:, self.n + 3] * np.vstack(aa[:, self.n + 3])

1187 ryy = aa[:, self.n + 4] * np.vstack(aa[:, self.n + 4])

1188 rzz = aa[:, self.n + 5] * np.vstack(aa[:, self.n + 5])

1189 T = txx + tyy + tzz + rxx + ryy + rzz

1190 for vec in hh:

1191 T += vec

1192 ff = np.dot(T, np.vstack(ff))

1193 forces[:, :] -= np.dot(T, np.vstack(ff)).reshape(-1, 3)

1194

1195 def __repr__(self):

1196 constraints = [repr(constr) for constr in self.constraints]

1197 return f'FixInternals(_copy_init={constraints}, epsilon={self.epsilon})'

1198

1199 # Classes for internal use in FixInternals

1200 class FixInternalsBase:

1201 """Base class for subclasses of FixInternals."""

1202

1203 def __init__(self, targetvalue, indices, masses, cell, pbc):

1204 self.targetvalue = targetvalue # constant target value

1205 self.indices = [defin[0:-1] for defin in indices] # indices, defs

1206 self.coefs = np.asarray([defin[-1] for defin in indices])

1207 self.masses = masses

1208 self.jacobian = [] # geometric Jacobian matrix, Wilson B-matrix

1209 self.sigma = 1. # difference between current and target value

1210 self.projected_force = None # helps optimizers scan along constr.

1211 self.cell = cell

1212 self.pbc = pbc

1213

1214 def finalize_jacobian(self, pos, n_internals, n, derivs):

1215 """Populate jacobian with derivatives for `n_internals` defined

1216 internals. n = 2 (bonds), 3 (angles), 4 (dihedrals)."""

1217 jacobian = np.zeros((n_internals, *pos.shape))

1218 for i, idx in enumerate(self.indices):

1219 for j in range(n):

1220 jacobian[i, idx[j]] = derivs[i, j]

1221 jacobian = jacobian.reshape((n_internals, 3 * len(pos)))

1222 return self.coefs @ jacobian

1223

1224 def finalize_positions(self, newpos):

1225 jacobian = self.jacobian / self.masses

1226 lamda = -self.sigma / (jacobian @ self.get_jacobian(newpos))

1227 dnewpos = lamda * jacobian

1228 newpos += dnewpos.reshape(newpos.shape)

1229

1230 def adjust_forces(self, positions, forces):

1231 self.projected_forces = ((self.jacobian @ forces.ravel())

1232 * self.jacobian)

1233 self.jacobian /= np.linalg.norm(self.jacobian)

1234

1235 class FixBondCombo(FixInternalsBase):

1236 """Constraint subobject for fixing linear combination of bond lengths

1237 within FixInternals.

1238

1239 sum_i( coef_i * bond_length_i ) = constant

1240 """

1241

1242 def get_jacobian(self, pos):

1243 bondvectors = [pos[k] - pos[h] for h, k in self.indices]

1244 derivs = get_distances_derivatives(bondvectors, cell=self.cell,

1245 pbc=self.pbc)

1246 return self.finalize_jacobian(pos, len(bondvectors), 2, derivs)

1247

1248 def setup_jacobian(self, pos):

1249 self.jacobian = self.get_jacobian(pos)

1250

1251 def adjust_positions(self, oldpos, newpos):

1252 bondvectors = [newpos[k] - newpos[h] for h, k in self.indices]

1253 (_, ), (dists, ) = conditional_find_mic([bondvectors],

1254 cell=self.cell,

1255 pbc=self.pbc)

1256 value = self.coefs @ dists

1257 self.sigma = value - self.targetvalue

1258 self.finalize_positions(newpos)

1259

1260 @staticmethod

1261 def get_value(atoms, indices, mic):

1262 return FixInternals.get_bondcombo(atoms, indices, mic)

1263

1264 def __repr__(self):

1265 return (f'FixBondCombo({self.targetvalue}, {self.indices}, '

1266 '{self.coefs})')

1267

1268 class FixBondLengthAlt(FixBondCombo):

1269 """Constraint subobject for fixing bond length within FixInternals.

1270 Fix distance between atoms with indices a1, a2."""

1271

1272 def __init__(self, targetvalue, indices, masses, cell, pbc):

1273 if targetvalue <= 0.:

1274 raise ZeroDivisionError('Invalid targetvalue for fixed bond')

1275 indices = [list(indices) + [1.]] # bond definition with coef 1.

1276 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)

1277

1278 @staticmethod

1279 def get_value(atoms, indices, mic):

1280 return atoms.get_distance(*indices, mic=mic)

1281

1282 def __repr__(self):

1283 return f'FixBondLengthAlt({self.targetvalue}, {self.indices})'

1284

1285 class FixAngle(FixInternalsBase):

1286 """Constraint subobject for fixing an angle within FixInternals.

1287

1288 Convergence is potentially problematic for angles very close to

1289 0 or 180 degrees as there is a singularity in the Cartesian derivative.

1290 Fixing planar angles is therefore not supported at the moment.

1291 """

1292

1293 def __init__(self, targetvalue, indices, masses, cell, pbc):

1294 """Fix atom movement to construct a constant angle."""

1295 if targetvalue <= 0. or targetvalue >= 180.:

1296 raise ZeroDivisionError('Invalid targetvalue for fixed angle')

1297 indices = [list(indices) + [1.]] # angle definition with coef 1.

1298 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)

1299

1300 def gather_vectors(self, pos):

1301 v0 = [pos[h] - pos[k] for h, k, l in self.indices]

1302 v1 = [pos[l] - pos[k] for h, k, l in self.indices]

1303 return v0, v1

1304

1305 def get_jacobian(self, pos):

1306 v0, v1 = self.gather_vectors(pos)

1307 derivs = get_angles_derivatives(v0, v1, cell=self.cell,

1308 pbc=self.pbc)

1309 return self.finalize_jacobian(pos, len(v0), 3, derivs)

1310

1311 def setup_jacobian(self, pos):

1312 self.jacobian = self.get_jacobian(pos)

1313

1314 def adjust_positions(self, oldpos, newpos):

1315 v0, v1 = self.gather_vectors(newpos)

1316 value = get_angles(v0, v1, cell=self.cell, pbc=self.pbc)

1317 self.sigma = value - self.targetvalue

1318 self.finalize_positions(newpos)

1319

1320 @staticmethod

1321 def get_value(atoms, indices, mic):

1322 return atoms.get_angle(*indices, mic=mic)

1323

1324 def __repr__(self):

1325 return f'FixAngle({self.targetvalue}, {self.indices})'

1326

1327 class FixDihedral(FixInternalsBase):

1328 """Constraint subobject for fixing a dihedral angle within FixInternals.

1329

1330 A dihedral becomes undefined when at least one of the inner two angles

1331 becomes planar. Make sure to avoid this situation.

1332 """

1333

1334 def __init__(self, targetvalue, indices, masses, cell, pbc):

1335 indices = [list(indices) + [1.]] # dihedral def. with coef 1.

1336 super().__init__(targetvalue, indices, masses, cell=cell, pbc=pbc)

1337

1338 def gather_vectors(self, pos):

1339 v0 = [pos[k] - pos[h] for h, k, l, m in self.indices]

1340 v1 = [pos[l] - pos[k] for h, k, l, m in self.indices]

1341 v2 = [pos[m] - pos[l] for h, k, l, m in self.indices]

1342 return v0, v1, v2

1343

1344 def get_jacobian(self, pos):

1345 v0, v1, v2 = self.gather_vectors(pos)

1346 derivs = get_dihedrals_derivatives(v0, v1, v2, cell=self.cell,

1347 pbc=self.pbc)

1348 return self.finalize_jacobian(pos, len(v0), 4, derivs)

1349

1350 def setup_jacobian(self, pos):

1351 self.jacobian = self.get_jacobian(pos)

1352

1353 def adjust_positions(self, oldpos, newpos):

1354 v0, v1, v2 = self.gather_vectors(newpos)

1355 value = get_dihedrals(v0, v1, v2, cell=self.cell, pbc=self.pbc)

1356 # apply minimum dihedral difference 'convention': (diff <= 180)

1357 self.sigma = (value - self.targetvalue + 180) % 360 - 180

1358 self.finalize_positions(newpos)

1359

1360 @staticmethod

1361 def get_value(atoms, indices, mic):

1362 return atoms.get_dihedral(*indices, mic=mic)

1363

1364 def __repr__(self):

1365 return f'FixDihedral({self.targetvalue}, {self.indices})'

1366

1367

1368class FixParametricRelations(FixConstraint):

1369

1370 def __init__(

1371 self,

1372 indices,

1373 Jacobian,

1374 const_shift,

1375 params=None,

1376 eps=1e-12,

1377 use_cell=False,

1378 ):

1379 """Constrains the degrees of freedom to act in a reduced parameter

1380 space defined by the Jacobian

1381

1382 These constraints are based off the work in:

1383 https://arxiv.org/abs/1908.01610

1384

1385 The constraints linearly maps the full 3N degrees of freedom,

1386 where N is number of active lattice vectors/atoms onto a

1387 reduced subset of M free parameters, where M <= 3*N. The

1388 Jacobian matrix and constant shift vector map the full set of

1389 degrees of freedom onto the reduced parameter space.

1390

1391 Currently the constraint is set up to handle either atomic

1392 positions or lattice vectors at one time, but not both. To do

1393 both simply add a two constraints for each set. This is done

1394 to keep the mathematics behind the operations separate.

1395

1396 It would be possible to extend these constraints to allow

1397 non-linear transformations if functionality to update the

1398 Jacobian at each position update was included. This would

1399 require passing an update function evaluate it every time

1400 adjust_positions is callled. This is currently NOT supported,

1401 and there are no plans to implement it in the future.

1402

1403 Args:

1404 indices (list of int): indices of the constrained atoms

1405 (if not None or empty then cell_indices must be None or Empty)

1406 Jacobian (np.ndarray(shape=(3*len(indices), len(params)))):

1407 The Jacobian describing

1408 the parameter space transformation

1409 const_shift (np.ndarray(shape=(3*len(indices)))):

1410 A vector describing the constant term

1411 in the transformation not accounted for in the Jacobian

1412 params (list of str):

1413 parameters used in the parametric representation

1414 if None a list is generated based on the shape of the Jacobian

1415 eps (float): a small number to compare the similarity of

1416 numbers and set the precision used

1417 to generate the constraint expressions

1418 use_cell (bool): if True then act on the cell object

1419

1420 """

1421 self.indices = np.array(indices)

1422 self.Jacobian = np.array(Jacobian)

1423 self.const_shift = np.array(const_shift)

1424

1425 assert self.const_shift.shape[0] == 3 * len(self.indices)

1426 assert self.Jacobian.shape[0] == 3 * len(self.indices)

1427

1428 self.eps = eps

1429 self.use_cell = use_cell

1430

1431 if params is None:

1432 params = []

1433 if self.Jacobian.shape[1] > 0:

1434 int_fmt_str = "{:0" + \

1435 str(int(np.ceil(np.log10(self.Jacobian.shape[1])))) + "d}"

1436 for param_ind in range(self.Jacobian.shape[1]):

1437 params.append("param_" + int_fmt_str.format(param_ind))

1438 else:

1439 assert len(params) == self.Jacobian.shape[-1]

1440

1441 self.params = params

1442

1443 self.Jacobian_inv = np.linalg.inv(

1444 self.Jacobian.T @ self.Jacobian) @ self.Jacobian.T

1445

1446 @classmethod

1447 def from_expressions(cls, indices, params, expressions,

1448 eps=1e-12, use_cell=False):

1449 """Converts the expressions into a Jacobian Matrix/const_shift

1450 vector and constructs a FixParametricRelations constraint

1451

1452 The expressions must be a list like object of size 3*N and

1453 elements must be ordered as:

1454 [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],

1455 where i, j, and k are the first, second and third

1456 component of the atomic position/lattice

1457 vector. Currently only linear operations are allowed to be

1458 included in the expressions so

1459 only terms like:

1460 - const * param_0

1461 - sqrt[const] * param_1

1462 - const * param_0 +/- const * param_1 +/- ... +/- const * param_M

1463 where const is any real number and param_0, param_1, ..., param_M are

1464 the parameters passed in

1465 params, are allowed.

1466

1467 For example, fractional atomic position constraints for wurtzite are:

1468 params = ["z1", "z2"]

1469 expressions = [

1470 "1.0/3.0", "2.0/3.0", "z1",

1471 "2.0/3.0", "1.0/3.0", "0.5 + z1",

1472 "1.0/3.0", "2.0/3.0", "z2",

1473 "2.0/3.0", "1.0/3.0", "0.5 + z2",

1474 ]

1475

1476 For diamond are:

1477 params = []

1478 expressions = [

1479 "0.0", "0.0", "0.0",

1480 "0.25", "0.25", "0.25",

1481 ],

1482

1483 and for stannite are

1484 params=["x4", "z4"]

1485 expressions = [

1486 "0.0", "0.0", "0.0",

1487 "0.0", "0.5", "0.5",

1488 "0.75", "0.25", "0.5",

1489 "0.25", "0.75", "0.5",

1490 "x4 + z4", "x4 + z4", "2*x4",

1491 "x4 - z4", "x4 - z4", "-2*x4",

1492 "0.0", "-1.0 * (x4 + z4)", "x4 - z4",

1493 "0.0", "x4 - z4", "-1.0 * (x4 + z4)",

1494 ]

1495

1496 Args:

1497 indices (list of int): indices of the constrained atoms

1498 (if not None or empty then cell_indices must be None or Empty)

1499 params (list of str): parameters used in the

1500 parametric representation

1501 expressions (list of str): expressions used to convert from the

1502 parametric to the real space representation

1503 eps (float): a small number to compare the similarity of

1504 numbers and set the precision used

1505 to generate the constraint expressions

1506 use_cell (bool): if True then act on the cell object

1507

1508 Returns:

1509 cls(

1510 indices,

1511 Jacobian generated from expressions,

1512 const_shift generated from expressions,

1513 params,

1514 eps-12,

1515 use_cell,

1516 )

1517 """

1518 Jacobian = np.zeros((3 * len(indices), len(params)))

1519 const_shift = np.zeros(3 * len(indices))

1520

1521 for expr_ind, expression in enumerate(expressions):

1522 expression = expression.strip()

1523

1524 # Convert subtraction to addition

1525 expression = expression.replace("-", "+(-1.0)*")

1526 if expression[0] == "+":

1527 expression = expression[1:]

1528 elif expression[:2] == "(+":

1529 expression = "(" + expression[2:]

1530

1531 # Explicitly add leading zeros so when replacing param_1 with 0.0

1532 # param_11 does not become 0.01

1533 int_fmt_str = "{:0" + \

1534 str(int(np.ceil(np.log10(len(params) + 1)))) + "d}"

1535

1536 param_dct = {}

1537 param_map = {}

1538

1539 # Construct a standardized param template for A/B filling

1540 for param_ind, param in enumerate(params):

1541 param_str = "param_" + int_fmt_str.format(param_ind)

1542 param_map[param] = param_str

1543 param_dct[param_str] = 0.0

1544

1545 # Replace the parameters according to the map

1546 # Sort by string length (long to short) to prevent cases like x11

1547 # becoming f"{param_map["x1"]}1"

1548 for param in sorted(params, key=lambda s: -1.0 * len(s)):

1549 expression = expression.replace(param, param_map[param])

1550

1551 # Partial linearity check

1552 for express_sec in expression.split("+"):

1553 in_sec = [param in express_sec for param in param_dct]

1554 n_params_in_sec = len(np.where(np.array(in_sec))[0])

1555 if n_params_in_sec > 1:

1556 raise ValueError(

1557 "FixParametricRelations expressions must be linear.")

1558

1559 const_shift[expr_ind] = float(

1560 eval_expression(expression, param_dct))

1561

1562 for param_ind in range(len(params)):

1563 param_str = "param_" + int_fmt_str.format(param_ind)

1564 if param_str not in expression:

1565 Jacobian[expr_ind, param_ind] = 0.0

1566 continue

1567 param_dct[param_str] = 1.0

1568 test_1 = float(eval_expression(expression, param_dct))

1569 test_1 -= const_shift[expr_ind]

1570 Jacobian[expr_ind, param_ind] = test_1

1571

1572 param_dct[param_str] = 2.0

1573 test_2 = float(eval_expression(expression, param_dct))

1574 test_2 -= const_shift[expr_ind]

1575 if abs(test_2 / test_1 - 2.0) > eps:

1576 raise ValueError(

1577 "FixParametricRelations expressions must be linear.")

1578 param_dct[param_str] = 0.0

1579

1580 args = [

1581 indices,

1582 Jacobian,

1583 const_shift,

1584 params,

1585 eps,

1586 use_cell,

1587 ]

1588 if cls is FixScaledParametricRelations:

1589 args = args[:-1]

1590 return cls(*args)

1591

1592 @property

1593 def expressions(self):

1594 """Generate the expressions represented by the current self.Jacobian

1595 and self.const_shift objects"""

1596 expressions = []

1597 per = int(round(-1 * np.log10(self.eps)))

1598 fmt_str = "{:." + str(per + 1) + "g}"

1599 for index, shift_val in enumerate(self.const_shift):

1600 exp = ""

1601 if np.all(np.abs(self.Jacobian[index]) < self.eps) or np.abs(

1602 shift_val) > self.eps:

1603 exp += fmt_str.format(shift_val)

1604

1605 param_exp = ""

1606 for param_index, jacob_val in enumerate(self.Jacobian[index]):

1607 abs_jacob_val = np.round(np.abs(jacob_val), per + 1)

1608 if abs_jacob_val < self.eps:

1609 continue

1610

1611 param = self.params[param_index]

1612 if param_exp or exp:

1613 if jacob_val > -1.0 * self.eps:

1614 param_exp += " + "

1615 else:

1616 param_exp += " - "

1617 elif (not exp) and (not param_exp) and (

1618 jacob_val < -1.0 * self.eps):

1619 param_exp += "-"

1620

1621 if np.abs(abs_jacob_val - 1.0) <= self.eps:

1622 param_exp += f"{param:s}"

1623 else:

1624 param_exp += (fmt_str +

1625 "*{:s}").format(abs_jacob_val, param)

1626

1627 exp += param_exp

1628

1629 expressions.append(exp)

1630 return np.array(expressions).reshape((-1, 3))

1631

1632 def todict(self):

1633 """Create a dictionary representation of the constraint"""

1634 return {

1635 "name": type(self).__name__,

1636 "kwargs": {

1637 "indices": self.indices,

1638 "params": self.params,

1639 "Jacobian": self.Jacobian,

1640 "const_shift": self.const_shift,

1641 "eps": self.eps,

1642 "use_cell": self.use_cell,

1643 }

1644 }

1645

1646 def __repr__(self):

1647 """The str representation of the constraint"""

1648 if len(self.indices) > 1:

1649 indices_str = "[{:d}, ..., {:d}]".format(

1650 self.indices[0], self.indices[-1])

1651 else:

1652 indices_str = f"[{self.indices[0]:d}]"

1653

1654 if len(self.params) > 1:

1655 params_str = "[{:s}, ..., {:s}]".format(

1656 self.params[0], self.params[-1])

1657 elif len(self.params) == 1:

1658 params_str = f"[{self.params[0]:s}]"

1659 else:

1660 params_str = "[]"

1661

1662 return '{:s}({:s}, {:s}, ..., {:e})'.format(

1663 type(self).__name__,

1664 indices_str,

1665 params_str,

1666 self.eps

1667 )

1668

1669

1670class FixScaledParametricRelations(FixParametricRelations):

1671

1672 def __init__(

1673 self,

1674 indices,

1675 Jacobian,

1676 const_shift,

1677 params=None,

1678 eps=1e-12,

1679 ):

1680 """The fractional coordinate version of FixParametricRelations

1681

1682 All arguments are the same, but since this is for fractional

1683 coordinates use_cell is false"""

1684 super().__init__(

1685 indices,

1686 Jacobian,

1687 const_shift,

1688 params,

1689 eps,

1690 False,

1691 )

1692

1693 def adjust_contravariant(self, cell, vecs, B):

1694 """Adjust the values of a set of vectors that are contravariant

1695 with the unit transformation"""

1696 scaled = cell.scaled_positions(vecs).flatten()

1697 scaled = self.Jacobian_inv @ (scaled - B)

1698 scaled = ((self.Jacobian @ scaled) + B).reshape((-1, 3))

1699

1700 return cell.cartesian_positions(scaled)

1701

1702 def adjust_positions(self, atoms, positions):

1703 """Adjust positions of the atoms to match the constraints"""

1704 positions[self.indices] = self.adjust_contravariant(

1705 atoms.cell,

1706 positions[self.indices],

1707 self.const_shift,

1708 )

1709 positions[self.indices] = self.adjust_B(

1710 atoms.cell, positions[self.indices])

1711

1712 def adjust_B(self, cell, positions):

1713 """Wraps the positions back to the unit cell and adjust B to

1714 keep track of this change"""

1715 fractional = cell.scaled_positions(positions)

1716 wrapped_fractional = (fractional % 1.0) % 1.0

1717 self.const_shift += np.round(wrapped_fractional - fractional).flatten()

1718 return cell.cartesian_positions(wrapped_fractional)

1719

1720 def adjust_momenta(self, atoms, momenta):

1721 """Adjust momenta of the atoms to match the constraints"""

1722 momenta[self.indices] = self.adjust_contravariant(

1723 atoms.cell,

1724 momenta[self.indices],

1725 np.zeros(self.const_shift.shape),

1726 )

1727

1728 def adjust_forces(self, atoms, forces):

1729 """Adjust forces of the atoms to match the constraints"""

1730 # Forces are coavarient to the coordinate transformation, use the

1731 # inverse transformations

1732 cart2frac_jacob = np.zeros(2 * (3 * len(atoms),))

1733 for i_atom in range(len(atoms)):

1734 cart2frac_jacob[3 * i_atom:3 * (i_atom + 1),

1735 3 * i_atom:3 * (i_atom + 1)] = atoms.cell.T

1736

1737 jacobian = cart2frac_jacob @ self.Jacobian

1738 jacobian_inv = np.linalg.inv(jacobian.T @ jacobian) @ jacobian.T

1739

1740 reduced_forces = jacobian.T @ forces.flatten()

1741 forces[self.indices] = (jacobian_inv.T @ reduced_forces).reshape(-1, 3)

1742

1743 def todict(self):

1744 """Create a dictionary representation of the constraint"""

1745 dct = super().todict()

1746 del dct["kwargs"]["use_cell"]

1747 return dct

1748

1749

1750class FixCartesianParametricRelations(FixParametricRelations):

1751

1752 def __init__(

1753 self,

1754 indices,

1755 Jacobian,

1756 const_shift,

1757 params=None,

1758 eps=1e-12,

1759 use_cell=False,

1760 ):

1761 """The Cartesian coordinate version of FixParametricRelations"""

1762 super().__init__(

1763 indices,

1764 Jacobian,

1765 const_shift,

1766 params,

1767 eps,

1768 use_cell,

1769 )

1770

1771 def adjust_contravariant(self, vecs, B):

1772 """Adjust the values of a set of vectors that are contravariant with

1773 the unit transformation"""

1774 vecs = self.Jacobian_inv @ (vecs.flatten() - B)

1775 vecs = ((self.Jacobian @ vecs) + B).reshape((-1, 3))

1776 return vecs

1777

1778 def adjust_positions(self, atoms, positions):

1779 """Adjust positions of the atoms to match the constraints"""

1780 if self.use_cell:

1781 return

1782 positions[self.indices] = self.adjust_contravariant(

1783 positions[self.indices],

1784 self.const_shift,

1785 )

1786

1787 def adjust_momenta(self, atoms, momenta):

1788 """Adjust momenta of the atoms to match the constraints"""

1789 if self.use_cell:

1790 return

1791 momenta[self.indices] = self.adjust_contravariant(

1792 momenta[self.indices],

1793 np.zeros(self.const_shift.shape),

1794 )

1795

1796 def adjust_forces(self, atoms, forces):

1797 """Adjust forces of the atoms to match the constraints"""

1798 if self.use_cell:

1799 return

1800

1801 forces_reduced = self.Jacobian.T @ forces[self.indices].flatten()

1802 forces[self.indices] = (self.Jacobian_inv.T @

1803 forces_reduced).reshape(-1, 3)

1804

1805 def adjust_cell(self, atoms, cell):

1806 """Adjust the cell of the atoms to match the constraints"""

1807 if not self.use_cell:

1808 return

1809 cell[self.indices] = self.adjust_contravariant(

1810 cell[self.indices],

1811 np.zeros(self.const_shift.shape),

1812 )

1813

1814 def adjust_stress(self, atoms, stress):

1815 """Adjust the stress of the atoms to match the constraints"""

1816 if not self.use_cell:

1817 return

1818

1819 stress_3x3 = voigt_6_to_full_3x3_stress(stress)

1820 stress_reduced = self.Jacobian.T @ stress_3x3[self.indices].flatten()

1821 stress_3x3[self.indices] = (

1822 self.Jacobian_inv.T @ stress_reduced).reshape(-1, 3)

1823

1824 stress[:] = full_3x3_to_voigt_6_stress(stress_3x3)

1825

1826

1827class Hookean(FixConstraint):

1828 """Applies a Hookean restorative force between a pair of atoms, an atom

1829 and a point, or an atom and a plane."""

1830

1831 def __init__(self, a1, a2, k, rt=None):

1832 """Forces two atoms to stay close together by applying no force if

1833 they are below a threshold length, rt, and applying a Hookean

1834 restorative force when the distance between them exceeds rt. Can

1835 also be used to tether an atom to a fixed point in space or to a

1836 distance above a plane.

1837

1838 a1 : int

1839 Index of atom 1

1840 a2 : one of three options

1841 1) index of atom 2

1842 2) a fixed point in cartesian space to which to tether a1

1843 3) a plane given as (A, B, C, D) in A x + B y + C z + D = 0.

1844 k : float

1845 Hooke's law (spring) constant to apply when distance

1846 exceeds threshold_length. Units of eV A^-2.

1847 rt : float

1848 The threshold length below which there is no force. The

1849 length is 1) between two atoms, 2) between atom and point.

1850 This argument is not supplied in case 3. Units of A.

1851

1852 If a plane is specified, the Hooke's law force is applied if the atom

1853 is on the normal side of the plane. For instance, the plane with

1854 (A, B, C, D) = (0, 0, 1, -7) defines a plane in the xy plane with a z

1855 intercept of +7 and a normal vector pointing in the +z direction.

1856 If the atom has z > 7, then a downward force would be applied of

1857 k * (atom.z - 7). The same plane with the normal vector pointing in

1858 the -z direction would be given by (A, B, C, D) = (0, 0, -1, 7).

1859

1860 References:

1861

1862 Andrew A. Peterson, Topics in Catalysis volume 57, pages40–53 (2014)

1863 https://link.springer.com/article/10.1007%2Fs11244-013-0161-8

1864 """

1865

1866 if isinstance(a2, int):

1867 self._type = 'two atoms'

1868 self.indices = [a1, a2]

1869 elif len(a2) == 3:

1870 self._type = 'point'

1871 self.index = a1

1872 self.origin = np.array(a2)

1873 elif len(a2) == 4:

1874 self._type = 'plane'

1875 self.index = a1

1876 self.plane = a2

1877 else:

1878 raise RuntimeError('Unknown type for a2')

1879 self.threshold = rt

1880 self.spring = k

1881

1882 def get_removed_dof(self, atoms):

1883 return 0

1884

1885 def todict(self):

1886 dct = {'name': 'Hookean'}

1887 dct['kwargs'] = {'rt': self.threshold,

1888 'k': self.spring}

1889 if self._type == 'two atoms':

1890 dct['kwargs']['a1'] = self.indices[0]

1891 dct['kwargs']['a2'] = self.indices[1]

1892 elif self._type == 'point':

1893 dct['kwargs']['a1'] = self.index

1894 dct['kwargs']['a2'] = self.origin

1895 elif self._type == 'plane':

1896 dct['kwargs']['a1'] = self.index

1897 dct['kwargs']['a2'] = self.plane

1898 else:

1899 raise NotImplementedError(f'Bad type: {self._type}')

1900 return dct

1901

1902 def adjust_positions(self, atoms, newpositions):

1903 pass

1904

1905 def adjust_momenta(self, atoms, momenta):

1906 pass

1907

1908 def adjust_forces(self, atoms, forces):

1909 positions = atoms.positions

1910 if self._type == 'plane':

1911 A, B, C, D = self.plane

1912 x, y, z = positions[self.index]

1913 d = ((A * x + B * y + C * z + D) /

1914 np.sqrt(A**2 + B**2 + C**2))

1915 if d < 0:

1916 return

1917 magnitude = self.spring * d

1918 direction = - np.array((A, B, C)) / np.linalg.norm((A, B, C))

1919 forces[self.index] += direction * magnitude

1920 return

1921 if self._type == 'two atoms':

1922 p1, p2 = positions[self.indices]

1923 elif self._type == 'point':

1924 p1 = positions[self.index]

1925 p2 = self.origin

1926 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)

1927 bondlength = np.linalg.norm(displace)

1928 if bondlength > self.threshold:

1929 magnitude = self.spring * (bondlength - self.threshold)

1930 direction = displace / np.linalg.norm(displace)

1931 if self._type == 'two atoms':

1932 forces[self.indices[0]] += direction * magnitude

1933 forces[self.indices[1]] -= direction * magnitude

1934 else:

1935 forces[self.index] += direction * magnitude

1936

1937 def adjust_potential_energy(self, atoms):

1938 """Returns the difference to the potential energy due to an active

1939 constraint. (That is, the quantity returned is to be added to the

1940 potential energy.)"""

1941 positions = atoms.positions

1942 if self._type == 'plane':

1943 A, B, C, D = self.plane

1944 x, y, z = positions[self.index]

1945 d = ((A * x + B * y + C * z + D) /

1946 np.sqrt(A**2 + B**2 + C**2))

1947 if d > 0:

1948 return 0.5 * self.spring * d**2

1949 else:

1950 return 0.

1951 if self._type == 'two atoms':

1952 p1, p2 = positions[self.indices]

1953 elif self._type == 'point':

1954 p1 = positions[self.index]

1955 p2 = self.origin

1956 displace, _ = find_mic(p2 - p1, atoms.cell, atoms.pbc)

1957 bondlength = np.linalg.norm(displace)

1958 if bondlength > self.threshold:

1959 return 0.5 * self.spring * (bondlength - self.threshold)**2

1960 else:

1961 return 0.

1962

1963 def get_indices(self):

1964 if self._type == 'two atoms':

1965 return self.indices

1966 elif self._type == 'point':

1967 return self.index

1968 elif self._type == 'plane':

1969 return self.index

1970

1971 def index_shuffle(self, atoms, ind):

1972 # See docstring of superclass

1973 if self._type == 'two atoms':

1974 newa = [-1, -1] # Signal error

1975 for new, old in slice2enlist(ind, len(atoms)):

1976 for i, a in enumerate(self.indices):

1977 if old == a:

1978 newa[i] = new

1979 if newa[0] == -1 or newa[1] == -1:

1980 raise IndexError('Constraint not part of slice')

1981 self.indices = newa

1982 elif (self._type == 'point') or (self._type == 'plane'):

1983 newa = -1 # Signal error

1984 for new, old in slice2enlist(ind, len(atoms)):

1985 if old == self.index:

1986 newa = new

1987 break

1988 if newa == -1:

1989 raise IndexError('Constraint not part of slice')

1990 self.index = newa

1991

1992 def __repr__(self):

1993 if self._type == 'two atoms':

1994 return 'Hookean(%d, %d)' % tuple(self.indices)

1995 elif self._type == 'point':

1996 return 'Hookean(%d) to cartesian' % self.index

1997 else:

1998 return 'Hookean(%d) to plane' % self.index

1999

2000

2001class ExternalForce(FixConstraint):

2002 """Constraint object for pulling two atoms apart by an external force.

2003

2004 You can combine this constraint for example with FixBondLength but make

2005 sure that *ExternalForce* comes first in the list if there are overlaps

2006 between atom1-2 and atom3-4:

2007

2008 >>> from ase.build import bulk

2009

2010 >>> atoms = bulk('Cu', 'fcc', a=3.6)

2011 >>> atom1, atom2, atom3, atom4 = atoms[:4]

2012 >>> fext = 1.0

2013 >>> con1 = ExternalForce(atom1, atom2, f_ext)

2014 >>> con2 = FixBondLength(atom3, atom4)

2015 >>> atoms.set_constraint([con1, con2])

2016

2017 see ase/test/external_force.py"""

2018

2019 def __init__(self, a1, a2, f_ext):

2020 self.indices = [a1, a2]

2021 self.external_force = f_ext

2022

2023 def get_removed_dof(self, atoms):

2024 return 0

2025

2026 def adjust_positions(self, atoms, new):

2027 pass

2028

2029 def adjust_forces(self, atoms, forces):

2030 dist = np.subtract.reduce(atoms.positions[self.indices])

2031 force = self.external_force * dist / np.linalg.norm(dist)

2032 forces[self.indices] += (force, -force)

2033

2034 def adjust_potential_energy(self, atoms):

2035 dist = np.subtract.reduce(atoms.positions[self.indices])

2036 return -np.linalg.norm(dist) * self.external_force

2037

2038 def index_shuffle(self, atoms, ind):

2039 """Shuffle the indices of the two atoms in this constraint"""

2040 newa = [-1, -1] # Signal error

2041 for new, old in slice2enlist(ind, len(atoms)):

2042 for i, a in enumerate(self.indices):

2043 if old == a:

2044 newa[i] = new

2045 if newa[0] == -1 or newa[1] == -1:

2046 raise IndexError('Constraint not part of slice')

2047 self.indices = newa

2048

2049 def __repr__(self):

2050 return 'ExternalForce(%d, %d, %f)' % (self.indices[0],

2051 self.indices[1],

2052 self.external_force)

2053

2054 def todict(self):

2055 return {'name': 'ExternalForce',

2056 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],

2057 'f_ext': self.external_force}}

2058

2059

2060class MirrorForce(FixConstraint):

2061 """Constraint object for mirroring the force between two atoms.

2062

2063 This class is designed to find a transition state with the help of a

2064 single optimization. It can be used if the transition state belongs to a

2065 bond breaking reaction. First the given bond length will be fixed until

2066 all other degrees of freedom are optimized, then the forces of the two

2067 atoms will be mirrored to find the transition state. The mirror plane is

2068 perpendicular to the connecting line of the atoms. Transition states in

2069 dependence of the force can be obtained by stretching the molecule and

2070 fixing its total length with *FixBondLength* or by using *ExternalForce*

2071 during the optimization with *MirrorForce*.

2072

2073 Parameters

2074 ----------

2075 a1: int

2076 First atom index.

2077 a2: int

2078 Second atom index.

2079 max_dist: float

2080 Upper limit of the bond length interval where the transition state

2081 can be found.

2082 min_dist: float

2083 Lower limit of the bond length interval where the transition state

2084 can be found.

2085 fmax: float

2086 Maximum force used for the optimization.

2087

2088 Notes

2089 -----

2090 You can combine this constraint for example with FixBondLength but make

2091 sure that *MirrorForce* comes first in the list if there are overlaps

2092 between atom1-2 and atom3-4:

2093

2094 >>> from ase.build import bulk

2095

2096 >>> atoms = bulk('Cu', 'fcc', a=3.6)

2097 >>> atom1, atom2, atom3, atom4 = atoms[:4]

2098 >>> con1 = MirrorForce(atom1, atom2)

2099 >>> con2 = FixBondLength(atom3, atom4)

2100 >>> atoms.set_constraint([con1, con2])

2101

2102 """

2103

2104 def __init__(self, a1, a2, max_dist=2.5, min_dist=1., fmax=0.1):

2105 self.indices = [a1, a2]

2106 self.min_dist = min_dist

2107 self.max_dist = max_dist

2108 self.fmax = fmax

2109

2110 def adjust_positions(self, atoms, new):

2111 pass

2112

2113 def adjust_forces(self, atoms, forces):

2114 dist = np.subtract.reduce(atoms.positions[self.indices])

2115 d = np.linalg.norm(dist)

2116 if (d < self.min_dist) or (d > self.max_dist):

2117 # Stop structure optimization

2118 forces[:] *= 0

2119 return

2120 dist /= d

2121 df = np.subtract.reduce(forces[self.indices])

2122 f = df.dot(dist)

2123 con_saved = atoms.constraints

2124 try:

2125 con = [con for con in con_saved

2126 if not isinstance(con, MirrorForce)]

2127 atoms.set_constraint(con)

2128 forces_copy = atoms.get_forces()

2129 finally:

2130 atoms.set_constraint(con_saved)

2131 df1 = -1 / 2. * f * dist

2132 forces_copy[self.indices] += (df1, -df1)

2133 # Check if forces would be converged if the bond with mirrored forces

2134 # would also be fixed

2135 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:

2136 factor = 1.

2137 else:

2138 factor = 0.

2139 df1 = -(1 + factor) / 2. * f * dist

2140 forces[self.indices] += (df1, -df1)

2141

2142 def index_shuffle(self, atoms, ind):

2143 """Shuffle the indices of the two atoms in this constraint

2144

2145 """

2146 newa = [-1, -1] # Signal error

2147 for new, old in slice2enlist(ind, len(atoms)):

2148 for i, a in enumerate(self.indices):

2149 if old == a:

2150 newa[i] = new

2151 if newa[0] == -1 or newa[1] == -1:

2152 raise IndexError('Constraint not part of slice')

2153 self.indices = newa

2154

2155 def __repr__(self):

2156 return 'MirrorForce(%d, %d, %f, %f, %f)' % (

2157 self.indices[0], self.indices[1], self.max_dist, self.min_dist,

2158 self.fmax)

2159

2160 def todict(self):

2161 return {'name': 'MirrorForce',

2162 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],

2163 'max_dist': self.max_dist,

2164 'min_dist': self.min_dist, 'fmax': self.fmax}}

2165

2166

2167class MirrorTorque(FixConstraint):

2168 """Constraint object for mirroring the torque acting on a dihedral

2169 angle defined by four atoms.

2170

2171 This class is designed to find a transition state with the help of a

2172 single optimization. It can be used if the transition state belongs to a

2173 cis-trans-isomerization with a change of dihedral angle. First the given

2174 dihedral angle will be fixed until all other degrees of freedom are

2175 optimized, then the torque acting on the dihedral angle will be mirrored

2176 to find the transition state. Transition states in

2177 dependence of the force can be obtained by stretching the molecule and

2178 fixing its total length with *FixBondLength* or by using *ExternalForce*

2179 during the optimization with *MirrorTorque*.

2180

2181 This constraint can be used to find

2182 transition states of cis-trans-isomerization.

2183

2184 a1 a4

2185 | |

2186 a2 __ a3

2187

2188 Parameters

2189 ----------

2190 a1: int

2191 First atom index.

2192 a2: int

2193 Second atom index.

2194 a3: int

2195 Third atom index.

2196 a4: int

2197 Fourth atom index.

2198 max_angle: float

2199 Upper limit of the dihedral angle interval where the transition state

2200 can be found.

2201 min_angle: float

2202 Lower limit of the dihedral angle interval where the transition state

2203 can be found.

2204 fmax: float

2205 Maximum force used for the optimization.

2206

2207 Notes

2208 -----

2209 You can combine this constraint for example with FixBondLength but make

2210 sure that *MirrorTorque* comes first in the list if there are overlaps

2211 between atom1-4 and atom5-6:

2212

2213 >>> from ase.build import bulk

2214

2215 >>> atoms = bulk('Cu', 'fcc', a=3.6)

2216 >>> atom1, atom2, atom3, atom4, atom5, atom6 = atoms[:6]

2217 >>> con1 = MirrorTorque(atom1, atom2, atom3, atom4)

2218 >>> con2 = FixBondLength(atom5, atom6)

2219 >>> atoms.set_constraint([con1, con2])

2220

2221 """

2222

2223 def __init__(self, a1, a2, a3, a4, max_angle=2 * np.pi, min_angle=0.,

2224 fmax=0.1):

2225 self.indices = [a1, a2, a3, a4]

2226 self.min_angle = min_angle

2227 self.max_angle = max_angle

2228 self.fmax = fmax

2229

2230 def adjust_positions(self, atoms, new):

2231 pass

2232

2233 def adjust_forces(self, atoms, forces):

2234 angle = atoms.get_dihedral(self.indices[0], self.indices[1],

2235 self.indices[2], self.indices[3])

2236 angle *= np.pi / 180.

2237 if (angle < self.min_angle) or (angle > self.max_angle):

2238 # Stop structure optimization

2239 forces[:] *= 0

2240 return

2241 p = atoms.positions[self.indices]

2242 f = forces[self.indices]

2243

2244 f0 = (f[1] + f[2]) / 2.

2245 ff = f - f0

2246 p0 = (p[2] + p[1]) / 2.

2247 m0 = np.cross(p[1] - p0, ff[1]) / (p[1] - p0).dot(p[1] - p0)

2248 fff = ff - np.cross(m0, p - p0)

2249 d1 = np.cross(np.cross(p[1] - p0, p[0] - p[1]), p[1] - p0) / \

2250 (p[1] - p0).dot(p[1] - p0)

2251 d2 = np.cross(np.cross(p[2] - p0, p[3] - p[2]), p[2] - p0) / \

2252 (p[2] - p0).dot(p[2] - p0)

2253 omegap1 = (np.cross(d1, fff[0]) / d1.dot(d1)).dot(p[1] - p0) / \

2254 np.linalg.norm(p[1] - p0)

2255 omegap2 = (np.cross(d2, fff[3]) / d2.dot(d2)).dot(p[2] - p0) / \

2256 np.linalg.norm(p[2] - p0)

2257 omegap = omegap1 + omegap2

2258 con_saved = atoms.constraints

2259 try:

2260 con = [con for con in con_saved

2261 if not isinstance(con, MirrorTorque)]

2262 atoms.set_constraint(con)

2263 forces_copy = atoms.get_forces()

2264 finally:

2265 atoms.set_constraint(con_saved)

2266 df1 = -1 / 2. * omegap * np.cross(p[1] - p0, d1) / \

2267 np.linalg.norm(p[1] - p0)

2268 df2 = -1 / 2. * omegap * np.cross(p[2] - p0, d2) / \

2269 np.linalg.norm(p[2] - p0)

2270 forces_copy[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)

2271 # Check if forces would be converged if the dihedral angle with

2272 # mirrored torque would also be fixed

2273 if (forces_copy**2).sum(axis=1).max() < self.fmax**2:

2274 factor = 1.

2275 else:

2276 factor = 0.

2277 df1 = -(1 + factor) / 2. * omegap * np.cross(p[1] - p0, d1) / \

2278 np.linalg.norm(p[1] - p0)

2279 df2 = -(1 + factor) / 2. * omegap * np.cross(p[2] - p0, d2) / \

2280 np.linalg.norm(p[2] - p0)

2281 forces[self.indices] += (df1, [0., 0., 0.], [0., 0., 0.], df2)

2282

2283 def index_shuffle(self, atoms, ind):

2284 # See docstring of superclass

2285 indices = []

2286 for new, old in slice2enlist(ind, len(atoms)):

2287 if old in self.indices:

2288 indices.append(new)

2289 if len(indices) == 0:

2290 raise IndexError('All indices in MirrorTorque not part of slice')

2291 self.indices = np.asarray(indices, int)

2292

2293 def __repr__(self):

2294 return 'MirrorTorque(%d, %d, %d, %d, %f, %f, %f)' % (

2295 self.indices[0], self.indices[1], self.indices[2],

2296 self.indices[3], self.max_angle, self.min_angle, self.fmax)

2297

2298 def todict(self):

2299 return {'name': 'MirrorTorque',

2300 'kwargs': {'a1': self.indices[0], 'a2': self.indices[1],

2301 'a3': self.indices[2], 'a4': self.indices[3],

2302 'max_angle': self.max_angle,

2303 'min_angle': self.min_angle, 'fmax': self.fmax}}

2304

2305

2306class FixSymmetry(FixConstraint):

2307 """

2308 Constraint to preserve spacegroup symmetry during optimisation.

2309

2310 Requires spglib package to be available.

2311 """

2312

2313 def __init__(self, atoms, symprec=0.01, adjust_positions=True,

2314 adjust_cell=True, verbose=False):

2315 self.atoms = atoms.copy()

2316 self.symprec = symprec

2317 self.verbose = verbose

2318 refine_symmetry(atoms, symprec, self.verbose) # refine initial symmetry

2319 sym = prep_symmetry(atoms, symprec, self.verbose)

2320 self.rotations, self.translations, self.symm_map = sym

2321 self.do_adjust_positions = adjust_positions

2322 self.do_adjust_cell = adjust_cell

2323

2324 def adjust_cell(self, atoms, cell):

2325 if not self.do_adjust_cell:

2326 return

2327 # stress should definitely be symmetrized as a rank 2 tensor

2328 # UnitCellFilter uses deformation gradient as cell DOF with steps

2329 # dF = stress.F^-T quantity that should be symmetrized is therefore dF .

2330 # F^T assume prev F = I, so just symmetrize dF

2331 cur_cell = atoms.get_cell()

2332 cur_cell_inv = atoms.cell.reciprocal().T

2333

2334 # F defined such that cell = cur_cell . F^T

2335 # assume prev F = I, so dF = F - I

2336 delta_deform_grad = np.dot(cur_cell_inv, cell).T - np.eye(3)

2337

2338 # symmetrization doesn't work properly with large steps, since

2339 # it depends on current cell, and cell is being changed by deformation

2340 # gradient

2341 max_delta_deform_grad = np.max(np.abs(delta_deform_grad))

2342 if max_delta_deform_grad > 0.25:

2343 raise RuntimeError('FixSymmetry adjust_cell does not work properly'

2344 ' with large deformation gradient step {} > 0.25'

2345 .format(max_delta_deform_grad))

2346 elif max_delta_deform_grad > 0.15:

2347 warn('FixSymmetry adjust_cell may be ill behaved with large '

2348 'deformation gradient step {}'.format(max_delta_deform_grad))

2349

2350 symmetrized_delta_deform_grad = symmetrize_rank2(cur_cell, cur_cell_inv,

2351 delta_deform_grad,

2352 self.rotations)

2353 cell[:] = np.dot(cur_cell,

2354 (symmetrized_delta_deform_grad + np.eye(3)).T)

2355

2356 def adjust_positions(self, atoms, new):

2357 if not self.do_adjust_positions:

2358 return

2359 # symmetrize changes in position as rank 1 tensors

2360 step = new - atoms.positions

2361 symmetrized_step = symmetrize_rank1(atoms.get_cell(),

2362 atoms.cell.reciprocal().T, step,

2363 self.rotations, self.translations,

2364 self.symm_map)

2365 new[:] = atoms.positions + symmetrized_step

2366

2367 def adjust_forces(self, atoms, forces):

2368 # symmetrize forces as rank 1 tensors

2369 # print('adjusting forces')

2370 forces[:] = symmetrize_rank1(atoms.get_cell(),

2371 atoms.cell.reciprocal().T, forces,

2372 self.rotations, self.translations,

2373 self.symm_map)

2374

2375 def adjust_stress(self, atoms, stress):

2376 # symmetrize stress as rank 2 tensor

2377 raw_stress = voigt_6_to_full_3x3_stress(stress)

2378 symmetrized_stress = symmetrize_rank2(atoms.get_cell(),

2379 atoms.cell.reciprocal().T,

2380 raw_stress, self.rotations)

2381 stress[:] = full_3x3_to_voigt_6_stress(symmetrized_stress)

2382

2383 def index_shuffle(self, atoms, ind):

2384 if len(atoms) != len(ind) or len(set(ind)) != len(ind):

2385 raise RuntimeError("FixSymmetry can only accomodate atom"

2386 " permutions, and len(Atoms) == {} "

2387 "!= len(ind) == {} or ind has duplicates"

2388 .format(len(atoms), len(ind)))

2389

2390 ind_reversed = np.zeros((len(ind)), dtype=int)

2391 ind_reversed[ind] = range(len(ind))

2392 new_symm_map = []

2393 for sm in self.symm_map:

2394 new_sm = np.array([-1] * len(atoms))

2395 for at_i in range(len(ind)):

2396 new_sm[ind_reversed[at_i]] = ind_reversed[sm[at_i]]

2397 new_symm_map.append(new_sm)

2398

2399 self.symm_map = new_symm_map

2400

2401 def todict(self):

2402 return {

2403 'name': 'FixSymmetry',

2404 'kwargs': {

2405 'atoms': self.atoms,

2406 'symprec': self.symprec,

2407 'adjust_positions': self.do_adjust_positions,

2408 'adjust_cell': self.do_adjust_cell,

2409 'verbose': self.verbose,

2410 },

2411 }

Coverage for ase / constraints.py: 87.36%

1203 statements