Coverage for /builds/kinetik161/ase/ase/filters.py: 86.24%
298 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Filters"""
2from itertools import product
3from warnings import warn
5import numpy as np
7from ase.calculators.calculator import PropertyNotImplementedError
8from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
9from ase.utils import deprecated, lazyproperty
10from ase.utils.abc import Optimizable
12__all__ = [
13 'Filter', 'StrainFilter', 'UnitCellFilter', 'FrechetCellFilter',
14 'ExpCellFilter'
15]
18class OptimizableFilter(Optimizable):
19 def __init__(self, filterobj):
20 self.filterobj = filterobj
22 def get_positions(self):
23 return self.filterobj.get_positions()
25 def set_positions(self, positions):
26 self.filterobj.set_positions(positions)
28 def get_forces(self):
29 return self.filterobj.get_forces()
31 @lazyproperty
32 def _use_force_consistent_energy(self):
33 # This boolean is in principle invalidated if the
34 # calculator changes. This can lead to weird things
35 # in multi-step optimizations.
36 try:
37 self.filterobj.get_potential_energy(force_consistent=True)
38 except PropertyNotImplementedError:
39 return False
40 else:
41 return True
43 def get_potential_energy(self):
44 force_consistent = self._use_force_consistent_energy
45 return self.filterobj.get_potential_energy(
46 force_consistent=force_consistent)
48 def __len__(self):
49 return len(self.filterobj)
51 def iterimages(self):
52 return self.filterobj.iterimages()
55class Filter:
56 """Subset filter class."""
58 def __init__(self, atoms, indices=None, mask=None):
59 """Filter atoms.
61 This filter can be used to hide degrees of freedom in an Atoms
62 object.
64 Parameters
65 ----------
66 indices : list of int
67 Indices for those atoms that should remain visible.
68 mask : list of bool
69 One boolean per atom indicating if the atom should remain
70 visible or not.
72 If a Trajectory tries to save this object, it will instead
73 save the underlying Atoms object. To prevent this, override
74 the iterimages method.
75 """
77 self.atoms = atoms
78 self.constraints = []
79 # Make self.info a reference to the underlying atoms' info dictionary.
80 self.info = self.atoms.info
82 if indices is None and mask is None:
83 raise ValueError('Use "indices" or "mask".')
84 if indices is not None and mask is not None:
85 raise ValueError('Use only one of "indices" and "mask".')
87 if mask is not None:
88 self.index = np.asarray(mask, bool)
89 self.n = self.index.sum()
90 else:
91 self.index = np.asarray(indices, int)
92 self.n = len(self.index)
94 def iterimages(self):
95 # Present the real atoms object to Trajectory and friends
96 return self.atoms.iterimages()
98 def get_cell(self):
99 """Returns the computational cell.
101 The computational cell is the same as for the original system.
102 """
103 return self.atoms.get_cell()
105 def get_pbc(self):
106 """Returns the periodic boundary conditions.
108 The boundary conditions are the same as for the original system.
109 """
110 return self.atoms.get_pbc()
112 def get_positions(self):
113 'Return the positions of the visible atoms.'
114 return self.atoms.get_positions()[self.index]
116 def set_positions(self, positions, **kwargs):
117 'Set the positions of the visible atoms.'
118 pos = self.atoms.get_positions()
119 pos[self.index] = positions
120 self.atoms.set_positions(pos, **kwargs)
122 positions = property(get_positions, set_positions,
123 doc='Positions of the atoms')
125 def get_momenta(self):
126 'Return the momenta of the visible atoms.'
127 return self.atoms.get_momenta()[self.index]
129 def set_momenta(self, momenta, **kwargs):
130 'Set the momenta of the visible atoms.'
131 mom = self.atoms.get_momenta()
132 mom[self.index] = momenta
133 self.atoms.set_momenta(mom, **kwargs)
135 def get_atomic_numbers(self):
136 'Return the atomic numbers of the visible atoms.'
137 return self.atoms.get_atomic_numbers()[self.index]
139 def set_atomic_numbers(self, atomic_numbers):
140 'Set the atomic numbers of the visible atoms.'
141 z = self.atoms.get_atomic_numbers()
142 z[self.index] = atomic_numbers
143 self.atoms.set_atomic_numbers(z)
145 def get_tags(self):
146 'Return the tags of the visible atoms.'
147 return self.atoms.get_tags()[self.index]
149 def set_tags(self, tags):
150 'Set the tags of the visible atoms.'
151 tg = self.atoms.get_tags()
152 tg[self.index] = tags
153 self.atoms.set_tags(tg)
155 def get_forces(self, *args, **kwargs):
156 return self.atoms.get_forces(*args, **kwargs)[self.index]
158 def get_stress(self, *args, **kwargs):
159 return self.atoms.get_stress(*args, **kwargs)
161 def get_stresses(self, *args, **kwargs):
162 return self.atoms.get_stresses(*args, **kwargs)[self.index]
164 def get_masses(self):
165 return self.atoms.get_masses()[self.index]
167 def get_potential_energy(self, **kwargs):
168 """Calculate potential energy.
170 Returns the potential energy of the full system.
171 """
172 return self.atoms.get_potential_energy(**kwargs)
174 def get_chemical_symbols(self):
175 return self.atoms.get_chemical_symbols()
177 def get_initial_magnetic_moments(self):
178 return self.atoms.get_initial_magnetic_moments()
180 def get_calculator(self):
181 """Returns the calculator.
183 WARNING: The calculator is unaware of this filter, and sees a
184 different number of atoms.
185 """
186 return self.atoms.calc
188 @property
189 def calc(self):
190 return self.atoms.calc
192 def get_celldisp(self):
193 return self.atoms.get_celldisp()
195 def has(self, name):
196 'Check for existence of array.'
197 return self.atoms.has(name)
199 def __len__(self):
200 'Return the number of movable atoms.'
201 return self.n
203 def __getitem__(self, i):
204 'Return an atom.'
205 return self.atoms[self.index[i]]
207 def __ase_optimizable__(self):
208 return OptimizableFilter(self)
211class StrainFilter(Filter):
212 """Modify the supercell while keeping the scaled positions fixed.
214 Presents the strain of the supercell as the generalized positions,
215 and the global stress tensor (times the volume) as the generalized
216 force.
218 This filter can be used to relax the unit cell until the stress is
219 zero. If MDMin is used for this, the timestep (dt) to be used
220 depends on the system size. 0.01/x where x is a typical dimension
221 seems like a good choice.
223 The stress and strain are presented as 6-vectors, the order of the
224 components follow the standard engingeering practice: xx, yy, zz,
225 yz, xz, xy.
227 """
229 def __init__(self, atoms, mask=None, include_ideal_gas=False):
230 """Create a filter applying a homogeneous strain to a list of atoms.
232 The first argument, atoms, is the atoms object.
234 The optional second argument, mask, is a list of six booleans,
235 indicating which of the six independent components of the
236 strain that are allowed to become non-zero. It defaults to
237 [1,1,1,1,1,1].
239 """
241 self.strain = np.zeros(6)
242 self.include_ideal_gas = include_ideal_gas
244 if mask is None:
245 mask = np.ones(6)
246 else:
247 mask = np.array(mask)
249 Filter.__init__(self, atoms=atoms, mask=mask)
250 self.mask = mask
251 self.origcell = atoms.get_cell()
253 def get_positions(self):
254 return self.strain.reshape((2, 3)).copy()
256 def set_positions(self, new):
257 new = new.ravel() * self.mask
258 eps = np.array([[1.0 + new[0], 0.5 * new[5], 0.5 * new[4]],
259 [0.5 * new[5], 1.0 + new[1], 0.5 * new[3]],
260 [0.5 * new[4], 0.5 * new[3], 1.0 + new[2]]])
262 self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True)
263 self.strain[:] = new
265 def get_forces(self, **kwargs):
266 stress = self.atoms.get_stress(include_ideal_gas=self.include_ideal_gas)
267 return -self.atoms.get_volume() * (stress * self.mask).reshape((2, 3))
269 def has(self, x):
270 return self.atoms.has(x)
272 def __len__(self):
273 return 2
276class UnitCellFilter(Filter):
277 """Modify the supercell and the atom positions. """
279 def __init__(self, atoms, mask=None,
280 cell_factor=None,
281 hydrostatic_strain=False,
282 constant_volume=False,
283 orig_cell=None,
284 scalar_pressure=0.0):
285 """Create a filter that returns the atomic forces and unit cell
286 stresses together, so they can simultaneously be minimized.
288 The first argument, atoms, is the atoms object. The optional second
289 argument, mask, is a list of booleans, indicating which of the six
290 independent components of the strain are relaxed.
292 - True = relax to zero
293 - False = fixed, ignore this component
295 Degrees of freedom are the positions in the original undeformed cell,
296 plus the deformation tensor (extra 3 "atoms"). This gives forces
297 consistent with numerical derivatives of the potential energy
298 with respect to the cell degreees of freedom.
300 For full details see:
301 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
302 Phys. Rev. B 59, 235 (1999)
304 You can still use constraints on the atoms, e.g. FixAtoms, to control
305 the relaxation of the atoms.
307 >>> # this should be equivalent to the StrainFilter
308 >>> atoms = Atoms(...)
309 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
310 >>> ucf = UnitCellFilter(atoms)
312 You should not attach this UnitCellFilter object to a
313 trajectory. Instead, create a trajectory for the atoms, and
314 attach it to an optimizer like this:
316 >>> atoms = Atoms(...)
317 >>> ucf = UnitCellFilter(atoms)
318 >>> qn = QuasiNewton(ucf)
319 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
320 >>> qn.attach(traj)
321 >>> qn.run(fmax=0.05)
323 Helpful conversion table:
325 - 0.05 eV/A^3 = 8 GPA
326 - 0.003 eV/A^3 = 0.48 GPa
327 - 0.0006 eV/A^3 = 0.096 GPa
328 - 0.0003 eV/A^3 = 0.048 GPa
329 - 0.0001 eV/A^3 = 0.02 GPa
331 Additional optional arguments:
333 cell_factor: float (default float(len(atoms)))
334 Factor by which deformation gradient is multiplied to put
335 it on the same scale as the positions when assembling
336 the combined position/cell vector. The stress contribution to
337 the forces is scaled down by the same factor. This can be thought
338 of as a very simple preconditioners. Default is number of atoms
339 which gives approximately the correct scaling.
341 hydrostatic_strain: bool (default False)
342 Constrain the cell by only allowing hydrostatic deformation.
343 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
345 constant_volume: bool (default False)
346 Project out the diagonal elements of the virial tensor to allow
347 relaxations at constant volume, e.g. for mapping out an
348 energy-volume curve. Note: this only approximately conserves
349 the volume and breaks energy/force consistency so can only be
350 used with optimizers that do require do a line minimisation
351 (e.g. FIRE).
353 scalar_pressure: float (default 0.0)
354 Applied pressure to use for enthalpy pV term. As above, this
355 breaks energy/force consistency.
356 """
358 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
359 self.atoms = atoms
360 if orig_cell is None:
361 self.orig_cell = atoms.get_cell()
362 else:
363 self.orig_cell = orig_cell
364 self.stress = None
366 if mask is None:
367 mask = np.ones(6)
368 mask = np.asarray(mask)
369 if mask.shape == (6,):
370 self.mask = voigt_6_to_full_3x3_stress(mask)
371 elif mask.shape == (3, 3):
372 self.mask = mask
373 else:
374 raise ValueError('shape of mask should be (3,3) or (6,)')
376 if cell_factor is None:
377 cell_factor = float(len(atoms))
378 self.hydrostatic_strain = hydrostatic_strain
379 self.constant_volume = constant_volume
380 self.scalar_pressure = scalar_pressure
381 self.cell_factor = cell_factor
382 self.copy = self.atoms.copy
383 self.arrays = self.atoms.arrays
385 def deform_grad(self):
386 return np.linalg.solve(self.orig_cell, self.atoms.cell).T
388 def get_positions(self):
389 """
390 this returns an array with shape (natoms + 3,3).
392 the first natoms rows are the positions of the atoms, the last
393 three rows are the deformation tensor associated with the unit cell,
394 scaled by self.cell_factor.
395 """
397 cur_deform_grad = self.deform_grad()
398 natoms = len(self.atoms)
399 pos = np.zeros((natoms + 3, 3))
400 # UnitCellFilter's positions are the self.atoms.positions but without
401 # the applied deformation gradient
402 pos[:natoms] = np.linalg.solve(cur_deform_grad,
403 self.atoms.positions.T).T
404 # UnitCellFilter's cell DOFs are the deformation gradient times a
405 # scaling factor
406 pos[natoms:] = self.cell_factor * cur_deform_grad
407 return pos
409 def set_positions(self, new, **kwargs):
410 """
411 new is an array with shape (natoms+3,3).
413 the first natoms rows are the positions of the atoms, the last
414 three rows are the deformation tensor used to change the cell shape.
416 the new cell is first set from original cell transformed by the new
417 deformation gradient, then the positions are set with respect to the
418 current cell by transforming them with the same deformation gradient
419 """
421 natoms = len(self.atoms)
422 new_atom_positions = new[:natoms]
423 new_deform_grad = new[natoms:] / self.cell_factor
424 # Set the new cell from the original cell and the new
425 # deformation gradient. Both current and final structures should
426 # preserve symmetry, so if set_cell() calls FixSymmetry.adjust_cell(),
427 # it should be OK
428 self.atoms.set_cell(self.orig_cell @ new_deform_grad.T,
429 scale_atoms=True)
430 # Set the positions from the ones passed in (which are without the
431 # deformation gradient applied) and the new deformation gradient.
432 # This should also preserve symmetry, so if set_positions() calls
433 # FixSymmetyr.adjust_positions(), it should be OK
434 self.atoms.set_positions(new_atom_positions @ new_deform_grad.T,
435 **kwargs)
437 def get_potential_energy(self, force_consistent=True):
438 """
439 returns potential energy including enthalpy PV term.
440 """
441 atoms_energy = self.atoms.get_potential_energy(
442 force_consistent=force_consistent)
443 return atoms_energy + self.scalar_pressure * self.atoms.get_volume()
445 def get_forces(self, **kwargs):
446 """
447 returns an array with shape (natoms+3,3) of the atomic forces
448 and unit cell stresses.
450 the first natoms rows are the forces on the atoms, the last
451 three rows are the forces on the unit cell, which are
452 computed from the stress tensor.
453 """
455 stress = self.atoms.get_stress(**kwargs)
456 atoms_forces = self.atoms.get_forces(**kwargs)
458 volume = self.atoms.get_volume()
459 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
460 np.diag([self.scalar_pressure] * 3))
461 cur_deform_grad = self.deform_grad()
462 atoms_forces = atoms_forces @ cur_deform_grad
463 virial = np.linalg.solve(cur_deform_grad, virial.T).T
465 if self.hydrostatic_strain:
466 vtr = virial.trace()
467 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
469 # Zero out components corresponding to fixed lattice elements
470 if (self.mask != 1.0).any():
471 virial *= self.mask
473 if self.constant_volume:
474 vtr = virial.trace()
475 np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0)
477 natoms = len(self.atoms)
478 forces = np.zeros((natoms + 3, 3))
479 forces[:natoms] = atoms_forces
480 forces[natoms:] = virial / self.cell_factor
482 self.stress = -full_3x3_to_voigt_6_stress(virial) / volume
483 return forces
485 def get_stress(self):
486 raise PropertyNotImplementedError
488 def has(self, x):
489 return self.atoms.has(x)
491 def __len__(self):
492 return (len(self.atoms) + 3)
495class FrechetCellFilter(UnitCellFilter):
496 """Modify the supercell and the atom positions."""
498 def __init__(self, atoms, mask=None,
499 exp_cell_factor=None,
500 hydrostatic_strain=False,
501 constant_volume=False,
502 scalar_pressure=0.0):
503 r"""Create a filter that returns the atomic forces and unit cell
504 stresses together, so they can simultaneously be minimized.
506 The first argument, atoms, is the atoms object. The optional second
507 argument, mask, is a list of booleans, indicating which of the six
508 independent components of the strain are relaxed.
510 - True = relax to zero
511 - False = fixed, ignore this component
513 Degrees of freedom are the positions in the original undeformed cell,
514 plus the log of the deformation tensor (extra 3 "atoms"). This gives
515 forces consistent with numerical derivatives of the potential energy
516 with respect to the cell degrees of freedom.
518 You can still use constraints on the atoms, e.g. FixAtoms, to control
519 the relaxation of the atoms.
521 >>> # this should be equivalent to the StrainFilter
522 >>> atoms = Atoms(...)
523 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
524 >>> ecf = FrechetCellFilter(atoms)
526 You should not attach this FrechetCellFilter object to a
527 trajectory. Instead, create a trajectory for the atoms, and
528 attach it to an optimizer like this:
530 >>> atoms = Atoms(...)
531 >>> ecf = FrechetCellFilter(atoms)
532 >>> qn = QuasiNewton(ecf)
533 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
534 >>> qn.attach(traj)
535 >>> qn.run(fmax=0.05)
537 Helpful conversion table:
539 - 0.05 eV/A^3 = 8 GPA
540 - 0.003 eV/A^3 = 0.48 GPa
541 - 0.0006 eV/A^3 = 0.096 GPa
542 - 0.0003 eV/A^3 = 0.048 GPa
543 - 0.0001 eV/A^3 = 0.02 GPa
545 Additional optional arguments:
547 exp_cell_factor: float (default float(len(atoms)))
548 Scaling factor for cell variables. The cell gradients in
549 FrechetCellFilter.get_forces() is divided by exp_cell_factor.
550 By default, set the number of atoms. We recommend to set
551 an extensive value for this parameter.
553 hydrostatic_strain: bool (default False)
554 Constrain the cell by only allowing hydrostatic deformation.
555 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
557 constant_volume: bool (default False)
558 Project out the diagonal elements of the virial tensor to allow
559 relaxations at constant volume, e.g. for mapping out an
560 energy-volume curve.
562 scalar_pressure: float (default 0.0)
563 Applied pressure to use for enthalpy pV term. As above, this
564 breaks energy/force consistency.
566 Implementation note:
568 The implementation is based on that of Christoph Ortner in JuLIP.jl:
569 https://github.com/JuliaMolSim/JuLIP.jl/blob/master/src/expcell.jl
571 The initial implementation of ExpCellFilter gave inconsistent gradients
572 for cell variables (matrix log of the deformation tensor). If you would
573 like to keep the previous behavior, please use ExpCellFilter.
575 The derivation of gradients of energy w.r.t positions and the log of the
576 deformation tensor is given in
577 https://github.com/lan496/lan496.github.io/blob/main/notes/cell_grad.pdf
578 """
580 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
581 UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
582 hydrostatic_strain=hydrostatic_strain,
583 constant_volume=constant_volume,
584 scalar_pressure=scalar_pressure)
586 # We defer the scipy import to avoid high immediate import overhead
587 from scipy.linalg import expm, expm_frechet, logm
588 self.expm = expm
589 self.logm = logm
590 self.expm_frechet = expm_frechet
592 # Scaling factor for cell gradients
593 if exp_cell_factor is None:
594 exp_cell_factor = float(len(atoms))
595 self.exp_cell_factor = exp_cell_factor
597 def get_positions(self):
598 pos = UnitCellFilter.get_positions(self)
599 natoms = len(self.atoms)
600 pos[natoms:] = self.logm(pos[natoms:]) * self.exp_cell_factor
601 return pos
603 def set_positions(self, new, **kwargs):
604 natoms = len(self.atoms)
605 new2 = new.copy()
606 new2[natoms:] = self.expm(new[natoms:] / self.exp_cell_factor)
607 UnitCellFilter.set_positions(self, new2, **kwargs)
609 def get_forces(self, **kwargs):
610 # forces on atoms are same as UnitCellFilter, we just
611 # need to modify the stress contribution
612 stress = self.atoms.get_stress(**kwargs)
613 volume = self.atoms.get_volume()
614 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
615 np.diag([self.scalar_pressure] * 3))
617 cur_deform_grad = self.deform_grad()
618 cur_deform_grad_log = self.logm(cur_deform_grad)
620 if self.hydrostatic_strain:
621 vtr = virial.trace()
622 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
624 # Zero out components corresponding to fixed lattice elements
625 if (self.mask != 1.0).any():
626 virial *= self.mask
628 # Cell gradient for UnitCellFilter
629 ucf_cell_grad = virial @ np.linalg.inv(cur_deform_grad.T)
631 # Cell gradient for FrechetCellFilter
632 deform_grad_log_force = np.zeros((3, 3))
633 for mu, nu in product(range(3), repeat=2):
634 dir = np.zeros((3, 3))
635 dir[mu, nu] = 1.0
636 # Directional derivative of deformation to (mu, nu) strain direction
637 expm_der = self.expm_frechet(
638 cur_deform_grad_log,
639 dir,
640 compute_expm=False
641 )
642 deform_grad_log_force[mu, nu] = np.sum(expm_der * ucf_cell_grad)
644 # Cauchy stress used for convergence testing
645 convergence_crit_stress = -(virial / volume)
646 if self.constant_volume:
647 # apply constraint to force
648 dglf_trace = deform_grad_log_force.trace()
649 np.fill_diagonal(deform_grad_log_force,
650 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
651 # apply constraint to Cauchy stress used for convergence testing
652 ccs_trace = convergence_crit_stress.trace()
653 np.fill_diagonal(convergence_crit_stress,
654 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
656 atoms_forces = self.atoms.get_forces(**kwargs)
657 atoms_forces = atoms_forces @ cur_deform_grad
659 # pack gradients into vector
660 natoms = len(self.atoms)
661 forces = np.zeros((natoms + 3, 3))
662 forces[:natoms] = atoms_forces
663 forces[natoms:] = deform_grad_log_force / self.exp_cell_factor
664 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
665 return forces
668class ExpCellFilter(UnitCellFilter):
670 @deprecated(DeprecationWarning(
671 'Use FrechetCellFilter for better convergence w.r.t. cell variables.'
672 ))
673 def __init__(self, atoms, mask=None,
674 cell_factor=None,
675 hydrostatic_strain=False,
676 constant_volume=False,
677 scalar_pressure=0.0):
678 r"""Create a filter that returns the atomic forces and unit cell
679 stresses together, so they can simultaneously be minimized.
681 The first argument, atoms, is the atoms object. The optional second
682 argument, mask, is a list of booleans, indicating which of the six
683 independent components of the strain are relaxed.
685 - True = relax to zero
686 - False = fixed, ignore this component
688 Degrees of freedom are the positions in the original undeformed
689 cell, plus the log of the deformation tensor (extra 3 "atoms"). This
690 gives forces consistent with numerical derivatives of the potential
691 energy with respect to the cell degrees of freedom.
693 For full details see:
694 E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
695 Phys. Rev. B 59, 235 (1999)
697 You can still use constraints on the atoms, e.g. FixAtoms, to
698 control the relaxation of the atoms.
700 >>> # this should be equivalent to the StrainFilter
701 >>> atoms = Atoms(...)
702 >>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
703 >>> ecf = ExpCellFilter(atoms)
705 You should not attach this ExpCellFilter object to a
706 trajectory. Instead, create a trajectory for the atoms, and
707 attach it to an optimizer like this:
709 >>> atoms = Atoms(...)
710 >>> ecf = ExpCellFilter(atoms)
711 >>> qn = QuasiNewton(ecf)
712 >>> traj = Trajectory('TiO2.traj', 'w', atoms)
713 >>> qn.attach(traj)
714 >>> qn.run(fmax=0.05)
716 Helpful conversion table:
718 - 0.05 eV/A^3 = 8 GPA
719 - 0.003 eV/A^3 = 0.48 GPa
720 - 0.0006 eV/A^3 = 0.096 GPa
721 - 0.0003 eV/A^3 = 0.048 GPa
722 - 0.0001 eV/A^3 = 0.02 GPa
724 Additional optional arguments:
726 cell_factor: (DEPRECATED)
727 Retained for backwards compatibility, but no longer used.
729 hydrostatic_strain: bool (default False)
730 Constrain the cell by only allowing hydrostatic deformation.
731 The virial tensor is replaced by np.diag([np.trace(virial)]*3).
733 constant_volume: bool (default False)
734 Project out the diagonal elements of the virial tensor to allow
735 relaxations at constant volume, e.g. for mapping out an
736 energy-volume curve.
738 scalar_pressure: float (default 0.0)
739 Applied pressure to use for enthalpy pV term. As above, this
740 breaks energy/force consistency.
742 Implementation details:
744 The implementation is based on that of Christoph Ortner in JuLIP.jl:
745 https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244
747 We decompose the deformation gradient as
749 F = exp(U) F0
750 x = F * F0^{-1} z = exp(U) z
752 If we write the energy as a function of U we can transform the
753 stress associated with a perturbation V into a derivative using a
754 linear map V -> L(U, V).
756 \phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) .
757 ( L(U, V) exp(-U) exp(U) z )
759 where
761 \nabla E(U) : V = [S exp(-U)'] : L(U,V)
762 = L'(U, S exp(-U)') : V
763 = L(U', S exp(-U)') : V
764 = L(U, S exp(-U)) : V (provided U = U')
766 where the : operator represents double contraction,
767 i.e. A:B = trace(A'B), and
769 F = deformation tensor - 3x3 matrix
770 F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here
771 U = cell degrees of freedom used here - 3x3 matrix
772 V = perturbation to cell DoFs - 3x3 matrix
773 v = perturbation to position DoFs
774 x = atomic positions in deformed cell
775 z = atomic positions in original cell
776 \phi = potential energy
777 S = stress tensor [3x3 matrix]
778 L(U, V) = directional derivative of exp at U in direction V, i.e
779 d/dt exp(U + t V)|_{t=0} = L(U, V)
781 This means we can write
783 d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V
785 and therefore the contribution to the gradient of the energy is
787 \nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij
789 .. deprecated:: 3.23.0
790 Use :class:`~ase.filters.FrechetCellFilter` for better convergence
791 w.r.t. cell variables.
792 """
793 Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
794 UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
795 cell_factor=cell_factor,
796 hydrostatic_strain=hydrostatic_strain,
797 constant_volume=constant_volume,
798 scalar_pressure=scalar_pressure)
799 if cell_factor is not None:
800 # cell_factor used in UnitCellFilter does not affect on gradients of
801 # ExpCellFilter.
802 warn("cell_factor is deprecated")
803 self.cell_factor = 1.0
805 # We defer the scipy import to avoid high immediate import overhead
806 from scipy.linalg import expm, logm
807 self.expm = expm
808 self.logm = logm
810 def get_forces(self, **kwargs):
811 forces = UnitCellFilter.get_forces(self, **kwargs)
813 # forces on atoms are same as UnitCellFilter, we just
814 # need to modify the stress contribution
815 stress = self.atoms.get_stress(**kwargs)
816 volume = self.atoms.get_volume()
817 virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
818 np.diag([self.scalar_pressure] * 3))
820 cur_deform_grad = self.deform_grad()
821 cur_deform_grad_log = self.logm(cur_deform_grad)
823 if self.hydrostatic_strain:
824 vtr = virial.trace()
825 virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
827 # Zero out components corresponding to fixed lattice elements
828 if (self.mask != 1.0).any():
829 virial *= self.mask
831 deform_grad_log_force_naive = virial.copy()
832 Y = np.zeros((6, 6))
833 Y[0:3, 0:3] = cur_deform_grad_log
834 Y[3:6, 3:6] = cur_deform_grad_log
835 Y[0:3, 3:6] = - virial @ self.expm(-cur_deform_grad_log)
836 deform_grad_log_force = -self.expm(Y)[0:3, 3:6]
837 for (i1, i2) in [(0, 1), (0, 2), (1, 2)]:
838 ff = 0.5 * (deform_grad_log_force[i1, i2] +
839 deform_grad_log_force[i2, i1])
840 deform_grad_log_force[i1, i2] = ff
841 deform_grad_log_force[i2, i1] = ff
843 # check for reasonable alignment between naive and
844 # exact search directions
845 all_are_equal = np.all(np.isclose(deform_grad_log_force,
846 deform_grad_log_force_naive))
847 if all_are_equal or \
848 (np.sum(deform_grad_log_force * deform_grad_log_force_naive) /
849 np.sqrt(np.sum(deform_grad_log_force**2) *
850 np.sum(deform_grad_log_force_naive**2)) > 0.8):
851 deform_grad_log_force = deform_grad_log_force_naive
853 # Cauchy stress used for convergence testing
854 convergence_crit_stress = -(virial / volume)
855 if self.constant_volume:
856 # apply constraint to force
857 dglf_trace = deform_grad_log_force.trace()
858 np.fill_diagonal(deform_grad_log_force,
859 np.diag(deform_grad_log_force) - dglf_trace / 3.0)
860 # apply constraint to Cauchy stress used for convergence testing
861 ccs_trace = convergence_crit_stress.trace()
862 np.fill_diagonal(convergence_crit_stress,
863 np.diag(convergence_crit_stress) - ccs_trace / 3.0)
865 # pack gradients into vector
866 natoms = len(self.atoms)
867 forces[natoms:] = deform_grad_log_force
868 self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
869 return forces