Coverage for /builds/ase/ase/ase/filters.py: 86.38%
301 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Filters"""
4from functools import cached_property
5from itertools import product
6from warnings import warn
8import numpy as np
10from ase.calculators.calculator import PropertyNotImplementedError
11from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
12from ase.utils import deprecated
13from ase.utils.abc import Optimizable
15__all__ = [
16 'Filter', 'StrainFilter', 'UnitCellFilter', 'FrechetCellFilter',
17 'ExpCellFilter'
18]
21class OptimizableFilter(Optimizable):
22 def __init__(self, filterobj):
23 self.filterobj = filterobj
25 def get_x(self):
26 return self.filterobj.get_positions().ravel()
28 def set_x(self, x):
29 self.filterobj.set_positions(x.reshape(-1, 3))
31 def get_gradient(self):
32 return self.filterobj.get_forces().ravel()
34 @cached_property
35 def _use_force_consistent_energy(self):
36 # This boolean is in principle invalidated if the
37 # calculator changes. This can lead to weird things
38 # in multi-step optimizations.
39 try:
40 self.filterobj.get_potential_energy(force_consistent=True)
41 except PropertyNotImplementedError:
42 return False
43 else:
44 return True
46 def get_value(self):
47 force_consistent = self._use_force_consistent_energy
48 return self.filterobj.get_potential_energy(
49 force_consistent=force_consistent)
51 def ndofs(self):
52 return 3 * len(self.filterobj)
54 def iterimages(self):
55 return self.filterobj.iterimages()
58class Filter:
59 """Subset filter class."""
61 def __init__(self, atoms, indices=None, mask=None):
62 """Filter atoms.
64 This filter can be used to hide degrees of freedom in an Atoms
65 object.
67 Parameters
68 ----------
69 indices : list of int
70 Indices for those atoms that should remain visible.
71 mask : list of bool
72 One boolean per atom indicating if the atom should remain
73 visible or not.
75 If a Trajectory tries to save this object, it will instead
76 save the underlying Atoms object. To prevent this, override
77 the iterimages method.
78 """
80 self.atoms = atoms
81 self.constraints = []
82 # Make self.info a reference to the underlying atoms' info dictionary.
83 self.info = self.atoms.info
85 if indices is None and mask is None:
86 raise ValueError('Use "indices" or "mask".')
87 if indices is not None and mask is not None:
88 raise ValueError('Use only one of "indices" and "mask".')
90 if mask is not None:
91 self.index = np.asarray(mask, bool)
92 self.n = self.index.sum()
93 else:
94 self.index = np.asarray(indices, int)
95 self.n = len(self.index)
97 def iterimages(self):
98 # Present the real atoms object to Trajectory and friends
99 return self.atoms.iterimages()
101 def get_cell(self):
102 """Returns the computational cell.
104 The computational cell is the same as for the original system.
105 """
106 return self.atoms.get_cell()
108 def get_pbc(self):
109 """Returns the periodic boundary conditions.
111 The boundary conditions are the same as for the original system.
112 """
113 return self.atoms.get_pbc()
115 def get_positions(self):
116 'Return the positions of the visible atoms.'
117 return self.atoms.get_positions()[self.index]
119 def set_positions(self, positions, **kwargs):
120 'Set the positions of the visible atoms.'
121 pos = self.atoms.get_positions()
122 pos[self.index] = positions
123 self.atoms.set_positions(pos, **kwargs)
125 positions = property(get_positions, set_positions,
126 doc='Positions of the atoms')
128 def get_momenta(self):
129 'Return the momenta of the visible atoms.'
130 return self.atoms.get_momenta()[self.index]
132 def set_momenta(self, momenta, **kwargs):
133 'Set the momenta of the visible atoms.'
134 mom = self.atoms.get_momenta()
135 mom[self.index] = momenta
136 self.atoms.set_momenta(mom, **kwargs)
138 def get_atomic_numbers(self):
139 'Return the atomic numbers of the visible atoms.'
140 return self.atoms.get_atomic_numbers()[self.index]
142 def set_atomic_numbers(self, atomic_numbers):
143 'Set the atomic numbers of the visible atoms.'
144 z = self.atoms.get_atomic_numbers()
145 z[self.index] = atomic_numbers
146 self.atoms.set_atomic_numbers(z)
148 def get_tags(self):
149 'Return the tags of the visible atoms.'
150 return self.atoms.get_tags()[self.index]
152 def set_tags(self, tags):
153 'Set the tags of the visible atoms.'
154 tg = self.atoms.get_tags()
155 tg[self.index] = tags
156 self.atoms.set_tags(tg)
158 def get_forces(self, *args, **kwargs):
159 return self.atoms.get_forces(*args, **kwargs)[self.index]
161 def get_stress(self, *args, **kwargs):
162 return self.atoms.get_stress(*args, **kwargs)
164 def get_stresses(self, *args, **kwargs):
165 return self.atoms.get_stresses(*args, **kwargs)[self.index]
167 def get_masses(self):
168 return self.atoms.get_masses()[self.index]
170 def get_potential_energy(self, **kwargs):
171 """Calculate potential energy.
173 Returns the potential energy of the full system.
174 """
175 return self.atoms.get_potential_energy(**kwargs)
177 def get_chemical_symbols(self):
178 return self.atoms.get_chemical_symbols()
180 def get_initial_magnetic_moments(self):
181 return self.atoms.get_initial_magnetic_moments()
183 def get_calculator(self):
184 """Returns the calculator.
186 WARNING: The calculator is unaware of this filter, and sees a
187 different number of atoms.
188 """
189 return self.atoms.calc
191 @property
192 def calc(self):
193 return self.atoms.calc
195 def get_celldisp(self):
196 return self.atoms.get_celldisp()
198 def has(self, name):
199 'Check for existence of array.'
200 return self.atoms.has(name)
202 def __len__(self):
203 'Return the number of movable atoms.'
204 return self.n
206 def __getitem__(self, i):
207 'Return an atom.'
208 return self.atoms[self.index[i]]
210 def __ase_optimizable__(self):
211 return OptimizableFilter(self)
214class StrainFilter(Filter):
215 """Modify the supercell while keeping the scaled positions fixed.
217 Presents the strain of the supercell as the generalized positions,
218 and the global stress tensor (times the volume) as the generalized
219 force.
221 This filter can be used to relax the unit cell until the stress is
222 zero. If MDMin is used for this, the timestep (dt) to be used
223 depends on the system size. 0.01/x where x is a typical dimension
224 seems like a good choice.
226 The stress and strain are presented as 6-vectors, the order of the
227 components follow the standard engingeering practice: xx, yy, zz,
228 yz, xz, xy.
230 """
232 def __init__(self, atoms, mask=None, include_ideal_gas=False):
233 """Create a filter applying a homogeneous strain to a list of atoms.
235 The first argument, atoms, is the atoms object.
237 The optional second argument, mask, is a list of six booleans,
238 indicating which of the six independent components of the
239 strain that are allowed to become non-zero. It defaults to
240 [1,1,1,1,1,1].
242 """
244 self.strain = np.zeros(6)
245 self.include_ideal_gas = include_ideal_gas
247 if mask is None:
248 mask = np.ones(6)
249 else:
250 mask = np.array(mask)
252 Filter.__init__(self, atoms=atoms, mask=mask)
253 self.mask = mask
254 self.origcell = atoms.get_cell()
256 def get_positions(self):
257 return self.strain.reshape((2, 3)).copy()
259 def set_positions(self, new):
260 new = new.ravel() * self.mask
261 eps = np.array([[1.0 + new[0], 0.5 * new[5], 0.5 * new[4]],
262 [0.5 * new[5], 1.0 + new[1], 0.5 * new[3]],
263 [0.5 * new[4], 0.5 * new[3], 1.0 + new[2]]])
265 self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True)
266 self.strain[:] = new
268 def get_forces(self, **kwargs):
269 stress = self.atoms.get_stress(include_ideal_gas=self.include_ideal_gas)
270 return -self.atoms.get_volume() * (stress * self.mask).reshape((2, 3))
272 def has(self, x):
273 return self.atoms.has(x)
275 def __len__(self):
276 return 2
279class UnitCellFilter(Filter):
280 """Modify the supercell and the atom positions. """
282 def __init__(self, atoms, mask=None,
283 cell_factor=None,
284 hydrostatic_strain=False,
285 constant_volume=False,
286 orig_cell=None,
287 scalar_pressure=0.0):
288 """Create a filter that returns the atomic forces and unit cell
289 stresses together, so they can simultaneously be minimized.
291 The first argument, atoms, is the atoms object. The optional second
292 argument, mask, is a list of booleans, indicating which of the six
293 independent components of the strain are relaxed.
295 - True = relax to zero
296 - False = fixed, ignore this component
298 Degrees of freedom are the positions in the original undeformed cell,
299 plus the deformation tensor (extra 3 "atoms"). This gives forces
300 consistent with numerical derivatives of the potential energy
301 with respect to the cell degreees of freedom.
303 For full details see:
304 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
305 Phys. Rev. B 59, 235 (1999)
307 You can still use constraints on the atoms, e.g. FixAtoms, to control
308 the relaxation of the atoms.
310 >>> # this should be equivalent to the StrainFilter
311 >>> atoms = Atoms(...)
312 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
313 >>> ucf = UnitCellFilter(atoms)
315 You should not attach this UnitCellFilter object to a
316 trajectory. Instead, create a trajectory for the atoms, and
317 attach it to an optimizer like this:
319 >>> atoms = Atoms(...)
320 >>> ucf = UnitCellFilter(atoms)
321 >>> qn = QuasiNewton(ucf)
322 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
323 >>> qn.attach(traj)
324 >>> qn.run(fmax=0.05)
326 Helpful conversion table:
328 - 0.05 eV/A^3 = 8 GPA
329 - 0.003 eV/A^3 = 0.48 GPa
330 - 0.0006 eV/A^3 = 0.096 GPa
331 - 0.0003 eV/A^3 = 0.048 GPa
332 - 0.0001 eV/A^3 = 0.02 GPa
334 Additional optional arguments:
336 cell_factor: float (default float(len(atoms)))
337 Factor by which deformation gradient is multiplied to put
338 it on the same scale as the positions when assembling
339 the combined position/cell vector. The stress contribution to
340 the forces is scaled down by the same factor. This can be thought
341 of as a very simple preconditioners. Default is number of atoms
342 which gives approximately the correct scaling.
344 hydrostatic_strain: bool (default False)
345 Constrain the cell by only allowing hydrostatic deformation.
346 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
348 constant_volume: bool (default False)
349 Project out the diagonal elements of the virial tensor to allow
350 relaxations at constant volume, e.g. for mapping out an
351 energy-volume curve. Note: this only approximately conserves
352 the volume and breaks energy/force consistency so can only be
353 used with optimizers that do require do a line minimisation
354 (e.g. FIRE).
356 scalar_pressure: float (default 0.0)
357 Applied pressure to use for enthalpy pV term. As above, this
358 breaks energy/force consistency.
359 """
361 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
362 self.atoms = atoms
363 if orig_cell is None:
364 self.orig_cell = atoms.get_cell()
365 else:
366 self.orig_cell = orig_cell
367 self.stress = None
369 if mask is None:
370 mask = np.ones(6)
371 mask = np.asarray(mask)
372 if mask.shape == (6,):
373 self.mask = voigt_6_to_full_3x3_stress(mask)
374 elif mask.shape == (3, 3):
375 self.mask = mask
376 else:
377 raise ValueError('shape of mask should be (3,3) or (6,)')
379 if cell_factor is None:
380 cell_factor = float(len(atoms))
381 self.hydrostatic_strain = hydrostatic_strain
382 self.constant_volume = constant_volume
383 self.scalar_pressure = scalar_pressure
384 self.cell_factor = cell_factor
385 self.copy = self.atoms.copy
386 self.arrays = self.atoms.arrays
388 def deform_grad(self):
389 return np.linalg.solve(self.orig_cell, self.atoms.cell).T
391 def get_positions(self):
392 """
393 this returns an array with shape (natoms + 3,3).
395 the first natoms rows are the positions of the atoms, the last
396 three rows are the deformation tensor associated with the unit cell,
397 scaled by self.cell_factor.
398 """
400 cur_deform_grad = self.deform_grad()
401 natoms = len(self.atoms)
402 pos = np.zeros((natoms + 3, 3))
403 # UnitCellFilter's positions are the self.atoms.positions but without
404 # the applied deformation gradient
405 pos[:natoms] = np.linalg.solve(cur_deform_grad,
406 self.atoms.positions.T).T
407 # UnitCellFilter's cell DOFs are the deformation gradient times a
408 # scaling factor
409 pos[natoms:] = self.cell_factor * cur_deform_grad
410 return pos
412 def set_positions(self, new, **kwargs):
413 """
414 new is an array with shape (natoms+3,3).
416 the first natoms rows are the positions of the atoms, the last
417 three rows are the deformation tensor used to change the cell shape.
419 the new cell is first set from original cell transformed by the new
420 deformation gradient, then the positions are set with respect to the
421 current cell by transforming them with the same deformation gradient
422 """
424 natoms = len(self.atoms)
425 new_atom_positions = new[:natoms]
426 new_deform_grad = new[natoms:] / self.cell_factor
427 deform = (new_deform_grad - np.eye(3)).T * self.mask
428 # Set the new cell from the original cell and the new
429 # deformation gradient. Both current and final structures should
430 # preserve symmetry, so if set_cell() calls FixSymmetry.adjust_cell(),
431 # it should be OK
432 newcell = self.orig_cell @ (np.eye(3) + deform)
434 self.atoms.set_cell(newcell,
435 scale_atoms=True)
436 # Set the positions from the ones passed in (which are without the
437 # deformation gradient applied) and the new deformation gradient.
438 # This should also preserve symmetry, so if set_positions() calls
439 # FixSymmetyr.adjust_positions(), it should be OK
440 self.atoms.set_positions(new_atom_positions @ (np.eye(3) + deform),
441 **kwargs)
443 def get_potential_energy(self, force_consistent=True):
444 """
445 returns potential energy including enthalpy PV term.
446 """
447 atoms_energy = self.atoms.get_potential_energy(
448 force_consistent=force_consistent)
449 return atoms_energy + self.scalar_pressure * self.atoms.get_volume()
451 def get_forces(self, **kwargs):
452 """
453 returns an array with shape (natoms+3,3) of the atomic forces
454 and unit cell stresses.
456 the first natoms rows are the forces on the atoms, the last
457 three rows are the forces on the unit cell, which are
458 computed from the stress tensor.
459 """
461 stress = self.atoms.get_stress(**kwargs)
462 atoms_forces = self.atoms.get_forces(**kwargs)
464 volume = self.atoms.get_volume()
465 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
466 np.diag([self.scalar_pressure] * 3))
467 cur_deform_grad = self.deform_grad()
468 atoms_forces = atoms_forces @ cur_deform_grad
469 virial = np.linalg.solve(cur_deform_grad, virial.T).T
471 if self.hydrostatic_strain:
472 vtr = virial.trace()
473 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
475 # Zero out components corresponding to fixed lattice elements
476 if (self.mask != 1.0).any():
477 virial *= self.mask
479 if self.constant_volume:
480 vtr = virial.trace()
481 np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0)
483 natoms = len(self.atoms)
484 forces = np.zeros((natoms + 3, 3))
485 forces[:natoms] = atoms_forces
486 forces[natoms:] = virial / self.cell_factor
488 self.stress = -full_3x3_to_voigt_6_stress(virial) / volume
489 return forces
491 def get_stress(self):
492 raise PropertyNotImplementedError
494 def has(self, x):
495 return self.atoms.has(x)
497 def __len__(self):
498 return (len(self.atoms) + 3)
501class FrechetCellFilter(UnitCellFilter):
502 """Modify the supercell and the atom positions."""
504 def __init__(self, atoms, mask=None,
505 exp_cell_factor=None,
506 hydrostatic_strain=False,
507 constant_volume=False,
508 scalar_pressure=0.0):
509 r"""Create a filter that returns the atomic forces and unit cell
510 stresses together, so they can simultaneously be minimized.
512 The first argument, atoms, is the atoms object. The optional second
513 argument, mask, is a list of booleans, indicating which of the six
514 independent components of the strain are relaxed.
516 - True = relax to zero
517 - False = fixed, ignore this component
519 Degrees of freedom are the positions in the original undeformed cell,
520 plus the log of the deformation tensor (extra 3 "atoms"). This gives
521 forces consistent with numerical derivatives of the potential energy
522 with respect to the cell degrees of freedom.
524 You can still use constraints on the atoms, e.g. FixAtoms, to control
525 the relaxation of the atoms.
527 >>> # this should be equivalent to the StrainFilter
528 >>> atoms = Atoms(...)
529 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
530 >>> ecf = FrechetCellFilter(atoms)
532 You should not attach this FrechetCellFilter object to a
533 trajectory. Instead, create a trajectory for the atoms, and
534 attach it to an optimizer like this:
536 >>> atoms = Atoms(...)
537 >>> ecf = FrechetCellFilter(atoms)
538 >>> qn = QuasiNewton(ecf)
539 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
540 >>> qn.attach(traj)
541 >>> qn.run(fmax=0.05)
543 Helpful conversion table:
545 - 0.05 eV/A^3 = 8 GPA
546 - 0.003 eV/A^3 = 0.48 GPa
547 - 0.0006 eV/A^3 = 0.096 GPa
548 - 0.0003 eV/A^3 = 0.048 GPa
549 - 0.0001 eV/A^3 = 0.02 GPa
551 Additional optional arguments:
553 exp_cell_factor: float (default float(len(atoms)))
554 Scaling factor for cell variables. The cell gradients in
555 FrechetCellFilter.get_forces() is divided by exp_cell_factor.
556 By default, set the number of atoms. We recommend to set
557 an extensive value for this parameter.
559 hydrostatic_strain: bool (default False)
560 Constrain the cell by only allowing hydrostatic deformation.
561 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
563 constant_volume: bool (default False)
564 Project out the diagonal elements of the virial tensor to allow
565 relaxations at constant volume, e.g. for mapping out an
566 energy-volume curve.
568 scalar_pressure: float (default 0.0)
569 Applied pressure to use for enthalpy pV term. As above, this
570 breaks energy/force consistency.
572 Implementation note:
574 The implementation is based on that of Christoph Ortner in JuLIP.jl:
575 https://github.com/JuliaMolSim/JuLIP.jl/blob/master/src/expcell.jl
577 The initial implementation of ExpCellFilter gave inconsistent gradients
578 for cell variables (matrix log of the deformation tensor). If you would
579 like to keep the previous behavior, please use ExpCellFilter.
581 The derivation of gradients of energy w.r.t positions and the log of the
582 deformation tensor is given in
583 https://github.com/lan496/lan496.github.io/blob/main/notes/cell_grad.pdf
584 """
586 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
587 UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
588 hydrostatic_strain=hydrostatic_strain,
589 constant_volume=constant_volume,
590 scalar_pressure=scalar_pressure)
592 # We defer the scipy import to avoid high immediate import overhead
593 from scipy.linalg import expm, expm_frechet, logm
594 self.expm = expm
595 self.logm = logm
596 self.expm_frechet = expm_frechet
598 # Scaling factor for cell gradients
599 if exp_cell_factor is None:
600 exp_cell_factor = float(len(atoms))
601 self.exp_cell_factor = exp_cell_factor
603 def get_positions(self):
604 pos = UnitCellFilter.get_positions(self)
605 natoms = len(self.atoms)
606 pos[natoms:] = self.logm(pos[natoms:]) * self.exp_cell_factor
607 return pos
609 def set_positions(self, new, **kwargs):
610 natoms = len(self.atoms)
611 new2 = new.copy()
612 new2[natoms:] = self.expm(new[natoms:] / self.exp_cell_factor)
613 UnitCellFilter.set_positions(self, new2, **kwargs)
615 def get_forces(self, **kwargs):
616 # forces on atoms are same as UnitCellFilter, we just
617 # need to modify the stress contribution
618 stress = self.atoms.get_stress(**kwargs)
619 volume = self.atoms.get_volume()
620 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
621 np.diag([self.scalar_pressure] * 3))
623 cur_deform_grad = self.deform_grad()
624 cur_deform_grad_log = self.logm(cur_deform_grad)
626 if self.hydrostatic_strain:
627 vtr = virial.trace()
628 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
630 # Zero out components corresponding to fixed lattice elements
631 if (self.mask != 1.0).any():
632 virial *= self.mask
634 # Cell gradient for UnitCellFilter
635 ucf_cell_grad = virial @ np.linalg.inv(cur_deform_grad.T)
637 # Cell gradient for FrechetCellFilter
638 deform_grad_log_force = np.zeros((3, 3))
639 for mu, nu in product(range(3), repeat=2):
640 dir = np.zeros((3, 3))
641 dir[mu, nu] = 1.0
642 # Directional derivative of deformation to (mu, nu) strain direction
643 expm_der = self.expm_frechet(
644 cur_deform_grad_log,
645 dir,
646 compute_expm=False
647 )
648 deform_grad_log_force[mu, nu] = np.sum(expm_der * ucf_cell_grad)
650 # Cauchy stress used for convergence testing
651 convergence_crit_stress = -(virial / volume)
652 if self.constant_volume:
653 # apply constraint to force
654 dglf_trace = deform_grad_log_force.trace()
655 np.fill_diagonal(deform_grad_log_force,
656 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
657 # apply constraint to Cauchy stress used for convergence testing
658 ccs_trace = convergence_crit_stress.trace()
659 np.fill_diagonal(convergence_crit_stress,
660 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
662 atoms_forces = self.atoms.get_forces(**kwargs)
663 atoms_forces = atoms_forces @ cur_deform_grad
665 # pack gradients into vector
666 natoms = len(self.atoms)
667 forces = np.zeros((natoms + 3, 3))
668 forces[:natoms] = atoms_forces
669 forces[natoms:] = deform_grad_log_force / self.exp_cell_factor
670 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
671 return forces
674class ExpCellFilter(UnitCellFilter):
676 @deprecated(DeprecationWarning(
677 'Use FrechetCellFilter for better convergence w.r.t. cell variables.'
678 ))
679 def __init__(self, atoms, mask=None,
680 cell_factor=None,
681 hydrostatic_strain=False,
682 constant_volume=False,
683 scalar_pressure=0.0):
684 r"""Create a filter that returns the atomic forces and unit cell
685 stresses together, so they can simultaneously be minimized.
687 The first argument, atoms, is the atoms object. The optional second
688 argument, mask, is a list of booleans, indicating which of the six
689 independent components of the strain are relaxed.
691 - True = relax to zero
692 - False = fixed, ignore this component
694 Degrees of freedom are the positions in the original undeformed
695 cell, plus the log of the deformation tensor (extra 3 "atoms"). This
696 gives forces consistent with numerical derivatives of the potential
697 energy with respect to the cell degrees of freedom.
699 For full details see:
700 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
701 Phys. Rev. B 59, 235 (1999)
703 You can still use constraints on the atoms, e.g. FixAtoms, to
704 control the relaxation of the atoms.
706 >>> # this should be equivalent to the StrainFilter
707 >>> atoms = Atoms(...)
708 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
709 >>> ecf = ExpCellFilter(atoms)
711 You should not attach this ExpCellFilter object to a
712 trajectory. Instead, create a trajectory for the atoms, and
713 attach it to an optimizer like this:
715 >>> atoms = Atoms(...)
716 >>> ecf = ExpCellFilter(atoms)
717 >>> qn = QuasiNewton(ecf)
718 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
719 >>> qn.attach(traj)
720 >>> qn.run(fmax=0.05)
722 Helpful conversion table:
724 - 0.05 eV/A^3 = 8 GPA
725 - 0.003 eV/A^3 = 0.48 GPa
726 - 0.0006 eV/A^3 = 0.096 GPa
727 - 0.0003 eV/A^3 = 0.048 GPa
728 - 0.0001 eV/A^3 = 0.02 GPa
730 Additional optional arguments:
732 cell_factor: (DEPRECATED)
733 Retained for backwards compatibility, but no longer used.
735 hydrostatic_strain: bool (default False)
736 Constrain the cell by only allowing hydrostatic deformation.
737 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
739 constant_volume: bool (default False)
740 Project out the diagonal elements of the virial tensor to allow
741 relaxations at constant volume, e.g. for mapping out an
742 energy-volume curve.
744 scalar_pressure: float (default 0.0)
745 Applied pressure to use for enthalpy pV term. As above, this
746 breaks energy/force consistency.
748 Implementation details:
750 The implementation is based on that of Christoph Ortner in JuLIP.jl:
751 https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244
753 We decompose the deformation gradient as
755 F = exp(U) F0
756 x = F * F0^{-1} z = exp(U) z
758 If we write the energy as a function of U we can transform the
759 stress associated with a perturbation V into a derivative using a
760 linear map V -> L(U, V).
762 \phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) .
763 ( L(U, V) exp(-U) exp(U) z )
765 where
767 \nabla E(U) : V = [S exp(-U)'] : L(U,V)
768 = L'(U, S exp(-U)') : V
769 = L(U', S exp(-U)') : V
770 = L(U, S exp(-U)) : V (provided U = U')
772 where the : operator represents double contraction,
773 i.e. A:B = trace(A'B), and
775 F = deformation tensor - 3x3 matrix
776 F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here
777 U = cell degrees of freedom used here - 3x3 matrix
778 V = perturbation to cell DoFs - 3x3 matrix
779 v = perturbation to position DoFs
780 x = atomic positions in deformed cell
781 z = atomic positions in original cell
782 \phi = potential energy
783 S = stress tensor [3x3 matrix]
784 L(U, V) = directional derivative of exp at U in direction V, i.e
785 d/dt exp(U + t V)|_{t=0} = L(U, V)
787 This means we can write
789 d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V
791 and therefore the contribution to the gradient of the energy is
793 \nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij
795 .. deprecated:: 3.23.0
796 Use :class:`~ase.filters.FrechetCellFilter` for better convergence
797 w.r.t. cell variables.
798 """
799 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
800 UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
801 cell_factor=cell_factor,
802 hydrostatic_strain=hydrostatic_strain,
803 constant_volume=constant_volume,
804 scalar_pressure=scalar_pressure)
805 if cell_factor is not None:
806 # cell_factor used in UnitCellFilter does not affect on gradients of
807 # ExpCellFilter.
808 warn("cell_factor is deprecated")
809 self.cell_factor = 1.0
811 # We defer the scipy import to avoid high immediate import overhead
812 from scipy.linalg import expm, logm
813 self.expm = expm
814 self.logm = logm
816 def get_forces(self, **kwargs):
817 forces = UnitCellFilter.get_forces(self, **kwargs)
819 # forces on atoms are same as UnitCellFilter, we just
820 # need to modify the stress contribution
821 stress = self.atoms.get_stress(**kwargs)
822 volume = self.atoms.get_volume()
823 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
824 np.diag([self.scalar_pressure] * 3))
826 cur_deform_grad = self.deform_grad()
827 cur_deform_grad_log = self.logm(cur_deform_grad)
829 if self.hydrostatic_strain:
830 vtr = virial.trace()
831 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
833 # Zero out components corresponding to fixed lattice elements
834 if (self.mask != 1.0).any():
835 virial *= self.mask
837 deform_grad_log_force_naive = virial.copy()
838 Y = np.zeros((6, 6))
839 Y[0:3, 0:3] = cur_deform_grad_log
840 Y[3:6, 3:6] = cur_deform_grad_log
841 Y[0:3, 3:6] = - virial @ self.expm(-cur_deform_grad_log)
842 deform_grad_log_force = -self.expm(Y)[0:3, 3:6]
843 for (i1, i2) in [(0, 1), (0, 2), (1, 2)]:
844 ff = 0.5 * (deform_grad_log_force[i1, i2] +
845 deform_grad_log_force[i2, i1])
846 deform_grad_log_force[i1, i2] = ff
847 deform_grad_log_force[i2, i1] = ff
849 # check for reasonable alignment between naive and
850 # exact search directions
851 all_are_equal = np.all(np.isclose(deform_grad_log_force,
852 deform_grad_log_force_naive))
853 if all_are_equal or \
854 (np.sum(deform_grad_log_force * deform_grad_log_force_naive) /
855 np.sqrt(np.sum(deform_grad_log_force**2) *
856 np.sum(deform_grad_log_force_naive**2)) > 0.8):
857 deform_grad_log_force = deform_grad_log_force_naive
859 # Cauchy stress used for convergence testing
860 convergence_crit_stress = -(virial / volume)
861 if self.constant_volume:
862 # apply constraint to force
863 dglf_trace = deform_grad_log_force.trace()
864 np.fill_diagonal(deform_grad_log_force,
865 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
866 # apply constraint to Cauchy stress used for convergence testing
867 ccs_trace = convergence_crit_stress.trace()
868 np.fill_diagonal(convergence_crit_stress,
869 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
871 # pack gradients into vector
872 natoms = len(self.atoms)
873 forces[natoms:] = deform_grad_log_force
874 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
875 return forces