"""Filters"""
from itertools import product
from warnings import warn
import numpy as np
from ase.calculators.calculator import PropertyNotImplementedError
from ase.stress import full_3x3_to_voigt_6_stress, voigt_6_to_full_3x3_stress
from ase.utils import deprecated, lazyproperty
from ase.utils.abc import Optimizable
__all__ = [
'Filter', 'StrainFilter', 'UnitCellFilter', 'FrechetCellFilter',
'ExpCellFilter'
]
class OptimizableFilter(Optimizable):
def __init__(self, filterobj):
self.filterobj = filterobj
def get_positions(self):
return self.filterobj.get_positions()
def set_positions(self, positions):
self.filterobj.set_positions(positions)
def get_forces(self):
return self.filterobj.get_forces()
@lazyproperty
def _use_force_consistent_energy(self):
# This boolean is in principle invalidated if the
# calculator changes. This can lead to weird things
# in multi-step optimizations.
try:
self.filterobj.get_potential_energy(force_consistent=True)
except PropertyNotImplementedError:
return False
else:
return True
def get_potential_energy(self):
force_consistent = self._use_force_consistent_energy
return self.filterobj.get_potential_energy(
force_consistent=force_consistent)
def __len__(self):
return len(self.filterobj)
def iterimages(self):
return self.filterobj.iterimages()
[docs]class Filter:
"""Subset filter class."""
def __init__(self, atoms, indices=None, mask=None):
"""Filter atoms.
This filter can be used to hide degrees of freedom in an Atoms
object.
Parameters
----------
indices : list of int
Indices for those atoms that should remain visible.
mask : list of bool
One boolean per atom indicating if the atom should remain
visible or not.
If a Trajectory tries to save this object, it will instead
save the underlying Atoms object. To prevent this, override
the iterimages method.
"""
self.atoms = atoms
self.constraints = []
# Make self.info a reference to the underlying atoms' info dictionary.
self.info = self.atoms.info
if indices is None and mask is None:
raise ValueError('Use "indices" or "mask".')
if indices is not None and mask is not None:
raise ValueError('Use only one of "indices" and "mask".')
if mask is not None:
self.index = np.asarray(mask, bool)
self.n = self.index.sum()
else:
self.index = np.asarray(indices, int)
self.n = len(self.index)
def iterimages(self):
# Present the real atoms object to Trajectory and friends
return self.atoms.iterimages()
def get_cell(self):
"""Returns the computational cell.
The computational cell is the same as for the original system.
"""
return self.atoms.get_cell()
def get_pbc(self):
"""Returns the periodic boundary conditions.
The boundary conditions are the same as for the original system.
"""
return self.atoms.get_pbc()
def get_positions(self):
'Return the positions of the visible atoms.'
return self.atoms.get_positions()[self.index]
def set_positions(self, positions, **kwargs):
'Set the positions of the visible atoms.'
pos = self.atoms.get_positions()
pos[self.index] = positions
self.atoms.set_positions(pos, **kwargs)
positions = property(get_positions, set_positions,
doc='Positions of the atoms')
def get_momenta(self):
'Return the momenta of the visible atoms.'
return self.atoms.get_momenta()[self.index]
def set_momenta(self, momenta, **kwargs):
'Set the momenta of the visible atoms.'
mom = self.atoms.get_momenta()
mom[self.index] = momenta
self.atoms.set_momenta(mom, **kwargs)
def get_atomic_numbers(self):
'Return the atomic numbers of the visible atoms.'
return self.atoms.get_atomic_numbers()[self.index]
def set_atomic_numbers(self, atomic_numbers):
'Set the atomic numbers of the visible atoms.'
z = self.atoms.get_atomic_numbers()
z[self.index] = atomic_numbers
self.atoms.set_atomic_numbers(z)
def get_tags(self):
'Return the tags of the visible atoms.'
return self.atoms.get_tags()[self.index]
def set_tags(self, tags):
'Set the tags of the visible atoms.'
tg = self.atoms.get_tags()
tg[self.index] = tags
self.atoms.set_tags(tg)
def get_forces(self, *args, **kwargs):
return self.atoms.get_forces(*args, **kwargs)[self.index]
def get_stress(self, *args, **kwargs):
return self.atoms.get_stress(*args, **kwargs)
def get_stresses(self, *args, **kwargs):
return self.atoms.get_stresses(*args, **kwargs)[self.index]
def get_masses(self):
return self.atoms.get_masses()[self.index]
def get_potential_energy(self, **kwargs):
"""Calculate potential energy.
Returns the potential energy of the full system.
"""
return self.atoms.get_potential_energy(**kwargs)
def get_chemical_symbols(self):
return self.atoms.get_chemical_symbols()
def get_initial_magnetic_moments(self):
return self.atoms.get_initial_magnetic_moments()
def get_calculator(self):
"""Returns the calculator.
WARNING: The calculator is unaware of this filter, and sees a
different number of atoms.
"""
return self.atoms.calc
@property
def calc(self):
return self.atoms.calc
def get_celldisp(self):
return self.atoms.get_celldisp()
def has(self, name):
'Check for existence of array.'
return self.atoms.has(name)
def __len__(self):
'Return the number of movable atoms.'
return self.n
def __getitem__(self, i):
'Return an atom.'
return self.atoms[self.index[i]]
def __ase_optimizable__(self):
return OptimizableFilter(self)
[docs]class StrainFilter(Filter):
"""Modify the supercell while keeping the scaled positions fixed.
Presents the strain of the supercell as the generalized positions,
and the global stress tensor (times the volume) as the generalized
force.
This filter can be used to relax the unit cell until the stress is
zero. If MDMin is used for this, the timestep (dt) to be used
depends on the system size. 0.01/x where x is a typical dimension
seems like a good choice.
The stress and strain are presented as 6-vectors, the order of the
components follow the standard engingeering practice: xx, yy, zz,
yz, xz, xy.
"""
def __init__(self, atoms, mask=None, include_ideal_gas=False):
"""Create a filter applying a homogeneous strain to a list of atoms.
The first argument, atoms, is the atoms object.
The optional second argument, mask, is a list of six booleans,
indicating which of the six independent components of the
strain that are allowed to become non-zero. It defaults to
[1,1,1,1,1,1].
"""
self.strain = np.zeros(6)
self.include_ideal_gas = include_ideal_gas
if mask is None:
mask = np.ones(6)
else:
mask = np.array(mask)
Filter.__init__(self, atoms=atoms, mask=mask)
self.mask = mask
self.origcell = atoms.get_cell()
def get_positions(self):
return self.strain.reshape((2, 3)).copy()
def set_positions(self, new):
new = new.ravel() * self.mask
eps = np.array([[1.0 + new[0], 0.5 * new[5], 0.5 * new[4]],
[0.5 * new[5], 1.0 + new[1], 0.5 * new[3]],
[0.5 * new[4], 0.5 * new[3], 1.0 + new[2]]])
self.atoms.set_cell(np.dot(self.origcell, eps), scale_atoms=True)
self.strain[:] = new
def get_forces(self, **kwargs):
stress = self.atoms.get_stress(include_ideal_gas=self.include_ideal_gas)
return -self.atoms.get_volume() * (stress * self.mask).reshape((2, 3))
def has(self, x):
return self.atoms.has(x)
def __len__(self):
return 2
[docs]class UnitCellFilter(Filter):
"""Modify the supercell and the atom positions. """
def __init__(self, atoms, mask=None,
cell_factor=None,
hydrostatic_strain=False,
constant_volume=False,
orig_cell=None,
scalar_pressure=0.0):
"""Create a filter that returns the atomic forces and unit cell
stresses together, so they can simultaneously be minimized.
The first argument, atoms, is the atoms object. The optional second
argument, mask, is a list of booleans, indicating which of the six
independent components of the strain are relaxed.
- True = relax to zero
- False = fixed, ignore this component
Degrees of freedom are the positions in the original undeformed cell,
plus the deformation tensor (extra 3 "atoms"). This gives forces
consistent with numerical derivatives of the potential energy
with respect to the cell degreees of freedom.
For full details see:
E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
Phys. Rev. B 59, 235 (1999)
You can still use constraints on the atoms, e.g. FixAtoms, to control
the relaxation of the atoms.
>>> # this should be equivalent to the StrainFilter
>>> atoms = Atoms(...)
>>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
>>> ucf = UnitCellFilter(atoms)
You should not attach this UnitCellFilter object to a
trajectory. Instead, create a trajectory for the atoms, and
attach it to an optimizer like this:
>>> atoms = Atoms(...)
>>> ucf = UnitCellFilter(atoms)
>>> qn = QuasiNewton(ucf)
>>> traj = Trajectory('TiO2.traj', 'w', atoms)
>>> qn.attach(traj)
>>> qn.run(fmax=0.05)
Helpful conversion table:
- 0.05 eV/A^3 = 8 GPA
- 0.003 eV/A^3 = 0.48 GPa
- 0.0006 eV/A^3 = 0.096 GPa
- 0.0003 eV/A^3 = 0.048 GPa
- 0.0001 eV/A^3 = 0.02 GPa
Additional optional arguments:
cell_factor: float (default float(len(atoms)))
Factor by which deformation gradient is multiplied to put
it on the same scale as the positions when assembling
the combined position/cell vector. The stress contribution to
the forces is scaled down by the same factor. This can be thought
of as a very simple preconditioners. Default is number of atoms
which gives approximately the correct scaling.
hydrostatic_strain: bool (default False)
Constrain the cell by only allowing hydrostatic deformation.
The virial tensor is replaced by np.diag([np.trace(virial)]*3).
constant_volume: bool (default False)
Project out the diagonal elements of the virial tensor to allow
relaxations at constant volume, e.g. for mapping out an
energy-volume curve. Note: this only approximately conserves
the volume and breaks energy/force consistency so can only be
used with optimizers that do require do a line minimisation
(e.g. FIRE).
scalar_pressure: float (default 0.0)
Applied pressure to use for enthalpy pV term. As above, this
breaks energy/force consistency.
"""
Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
self.atoms = atoms
if orig_cell is None:
self.orig_cell = atoms.get_cell()
else:
self.orig_cell = orig_cell
self.stress = None
if mask is None:
mask = np.ones(6)
mask = np.asarray(mask)
if mask.shape == (6,):
self.mask = voigt_6_to_full_3x3_stress(mask)
elif mask.shape == (3, 3):
self.mask = mask
else:
raise ValueError('shape of mask should be (3,3) or (6,)')
if cell_factor is None:
cell_factor = float(len(atoms))
self.hydrostatic_strain = hydrostatic_strain
self.constant_volume = constant_volume
self.scalar_pressure = scalar_pressure
self.cell_factor = cell_factor
self.copy = self.atoms.copy
self.arrays = self.atoms.arrays
def deform_grad(self):
return np.linalg.solve(self.orig_cell, self.atoms.cell).T
def get_positions(self):
"""
this returns an array with shape (natoms + 3,3).
the first natoms rows are the positions of the atoms, the last
three rows are the deformation tensor associated with the unit cell,
scaled by self.cell_factor.
"""
cur_deform_grad = self.deform_grad()
natoms = len(self.atoms)
pos = np.zeros((natoms + 3, 3))
# UnitCellFilter's positions are the self.atoms.positions but without
# the applied deformation gradient
pos[:natoms] = np.linalg.solve(cur_deform_grad,
self.atoms.positions.T).T
# UnitCellFilter's cell DOFs are the deformation gradient times a
# scaling factor
pos[natoms:] = self.cell_factor * cur_deform_grad
return pos
def set_positions(self, new, **kwargs):
"""
new is an array with shape (natoms+3,3).
the first natoms rows are the positions of the atoms, the last
three rows are the deformation tensor used to change the cell shape.
the new cell is first set from original cell transformed by the new
deformation gradient, then the positions are set with respect to the
current cell by transforming them with the same deformation gradient
"""
natoms = len(self.atoms)
new_atom_positions = new[:natoms]
new_deform_grad = new[natoms:] / self.cell_factor
# Set the new cell from the original cell and the new
# deformation gradient. Both current and final structures should
# preserve symmetry, so if set_cell() calls FixSymmetry.adjust_cell(),
# it should be OK
self.atoms.set_cell(self.orig_cell @ new_deform_grad.T,
scale_atoms=True)
# Set the positions from the ones passed in (which are without the
# deformation gradient applied) and the new deformation gradient.
# This should also preserve symmetry, so if set_positions() calls
# FixSymmetyr.adjust_positions(), it should be OK
self.atoms.set_positions(new_atom_positions @ new_deform_grad.T,
**kwargs)
def get_potential_energy(self, force_consistent=True):
"""
returns potential energy including enthalpy PV term.
"""
atoms_energy = self.atoms.get_potential_energy(
force_consistent=force_consistent)
return atoms_energy + self.scalar_pressure * self.atoms.get_volume()
def get_forces(self, **kwargs):
"""
returns an array with shape (natoms+3,3) of the atomic forces
and unit cell stresses.
the first natoms rows are the forces on the atoms, the last
three rows are the forces on the unit cell, which are
computed from the stress tensor.
"""
stress = self.atoms.get_stress(**kwargs)
atoms_forces = self.atoms.get_forces(**kwargs)
volume = self.atoms.get_volume()
virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
np.diag([self.scalar_pressure] * 3))
cur_deform_grad = self.deform_grad()
atoms_forces = atoms_forces @ cur_deform_grad
virial = np.linalg.solve(cur_deform_grad, virial.T).T
if self.hydrostatic_strain:
vtr = virial.trace()
virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
# Zero out components corresponding to fixed lattice elements
if (self.mask != 1.0).any():
virial *= self.mask
if self.constant_volume:
vtr = virial.trace()
np.fill_diagonal(virial, np.diag(virial) - vtr / 3.0)
natoms = len(self.atoms)
forces = np.zeros((natoms + 3, 3))
forces[:natoms] = atoms_forces
forces[natoms:] = virial / self.cell_factor
self.stress = -full_3x3_to_voigt_6_stress(virial) / volume
return forces
def get_stress(self):
raise PropertyNotImplementedError
def has(self, x):
return self.atoms.has(x)
def __len__(self):
return (len(self.atoms) + 3)
[docs]class FrechetCellFilter(UnitCellFilter):
"""Modify the supercell and the atom positions."""
def __init__(self, atoms, mask=None,
exp_cell_factor=None,
hydrostatic_strain=False,
constant_volume=False,
scalar_pressure=0.0):
r"""Create a filter that returns the atomic forces and unit cell
stresses together, so they can simultaneously be minimized.
The first argument, atoms, is the atoms object. The optional second
argument, mask, is a list of booleans, indicating which of the six
independent components of the strain are relaxed.
- True = relax to zero
- False = fixed, ignore this component
Degrees of freedom are the positions in the original undeformed cell,
plus the log of the deformation tensor (extra 3 "atoms"). This gives
forces consistent with numerical derivatives of the potential energy
with respect to the cell degrees of freedom.
You can still use constraints on the atoms, e.g. FixAtoms, to control
the relaxation of the atoms.
>>> # this should be equivalent to the StrainFilter
>>> atoms = Atoms(...)
>>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
>>> ecf = FrechetCellFilter(atoms)
You should not attach this FrechetCellFilter object to a
trajectory. Instead, create a trajectory for the atoms, and
attach it to an optimizer like this:
>>> atoms = Atoms(...)
>>> ecf = FrechetCellFilter(atoms)
>>> qn = QuasiNewton(ecf)
>>> traj = Trajectory('TiO2.traj', 'w', atoms)
>>> qn.attach(traj)
>>> qn.run(fmax=0.05)
Helpful conversion table:
- 0.05 eV/A^3 = 8 GPA
- 0.003 eV/A^3 = 0.48 GPa
- 0.0006 eV/A^3 = 0.096 GPa
- 0.0003 eV/A^3 = 0.048 GPa
- 0.0001 eV/A^3 = 0.02 GPa
Additional optional arguments:
exp_cell_factor: float (default float(len(atoms)))
Scaling factor for cell variables. The cell gradients in
FrechetCellFilter.get_forces() is divided by exp_cell_factor.
By default, set the number of atoms. We recommend to set
an extensive value for this parameter.
hydrostatic_strain: bool (default False)
Constrain the cell by only allowing hydrostatic deformation.
The virial tensor is replaced by np.diag([np.trace(virial)]*3).
constant_volume: bool (default False)
Project out the diagonal elements of the virial tensor to allow
relaxations at constant volume, e.g. for mapping out an
energy-volume curve.
scalar_pressure: float (default 0.0)
Applied pressure to use for enthalpy pV term. As above, this
breaks energy/force consistency.
Implementation note:
The implementation is based on that of Christoph Ortner in JuLIP.jl:
https://github.com/JuliaMolSim/JuLIP.jl/blob/master/src/expcell.jl
The initial implementation of ExpCellFilter gave inconsistent gradients
for cell variables (matrix log of the deformation tensor). If you would
like to keep the previous behavior, please use ExpCellFilter.
The derivation of gradients of energy w.r.t positions and the log of the
deformation tensor is given in
https://github.com/lan496/lan496.github.io/blob/main/notes/cell_grad.pdf
"""
Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
hydrostatic_strain=hydrostatic_strain,
constant_volume=constant_volume,
scalar_pressure=scalar_pressure)
# We defer the scipy import to avoid high immediate import overhead
from scipy.linalg import expm, expm_frechet, logm
self.expm = expm
self.logm = logm
self.expm_frechet = expm_frechet
# Scaling factor for cell gradients
if exp_cell_factor is None:
exp_cell_factor = float(len(atoms))
self.exp_cell_factor = exp_cell_factor
def get_positions(self):
pos = UnitCellFilter.get_positions(self)
natoms = len(self.atoms)
pos[natoms:] = self.logm(pos[natoms:]) * self.exp_cell_factor
return pos
def set_positions(self, new, **kwargs):
natoms = len(self.atoms)
new2 = new.copy()
new2[natoms:] = self.expm(new[natoms:] / self.exp_cell_factor)
UnitCellFilter.set_positions(self, new2, **kwargs)
def get_forces(self, **kwargs):
# forces on atoms are same as UnitCellFilter, we just
# need to modify the stress contribution
stress = self.atoms.get_stress(**kwargs)
volume = self.atoms.get_volume()
virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
np.diag([self.scalar_pressure] * 3))
cur_deform_grad = self.deform_grad()
cur_deform_grad_log = self.logm(cur_deform_grad)
if self.hydrostatic_strain:
vtr = virial.trace()
virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
# Zero out components corresponding to fixed lattice elements
if (self.mask != 1.0).any():
virial *= self.mask
# Cell gradient for UnitCellFilter
ucf_cell_grad = virial @ np.linalg.inv(cur_deform_grad.T)
# Cell gradient for FrechetCellFilter
deform_grad_log_force = np.zeros((3, 3))
for mu, nu in product(range(3), repeat=2):
dir = np.zeros((3, 3))
dir[mu, nu] = 1.0
# Directional derivative of deformation to (mu, nu) strain direction
expm_der = self.expm_frechet(
cur_deform_grad_log,
dir,
compute_expm=False
)
deform_grad_log_force[mu, nu] = np.sum(expm_der * ucf_cell_grad)
# Cauchy stress used for convergence testing
convergence_crit_stress = -(virial / volume)
if self.constant_volume:
# apply constraint to force
dglf_trace = deform_grad_log_force.trace()
np.fill_diagonal(deform_grad_log_force,
np.diag(deform_grad_log_force) - dglf_trace / 3.0)
# apply constraint to Cauchy stress used for convergence testing
ccs_trace = convergence_crit_stress.trace()
np.fill_diagonal(convergence_crit_stress,
np.diag(convergence_crit_stress) - ccs_trace / 3.0)
atoms_forces = self.atoms.get_forces(**kwargs)
atoms_forces = atoms_forces @ cur_deform_grad
# pack gradients into vector
natoms = len(self.atoms)
forces = np.zeros((natoms + 3, 3))
forces[:natoms] = atoms_forces
forces[natoms:] = deform_grad_log_force / self.exp_cell_factor
self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
return forces
[docs]class ExpCellFilter(UnitCellFilter):
@deprecated(DeprecationWarning(
'Use FrechetCellFilter for better convergence w.r.t. cell variables.'
))
def __init__(self, atoms, mask=None,
cell_factor=None,
hydrostatic_strain=False,
constant_volume=False,
scalar_pressure=0.0):
r"""Create a filter that returns the atomic forces and unit cell
stresses together, so they can simultaneously be minimized.
The first argument, atoms, is the atoms object. The optional second
argument, mask, is a list of booleans, indicating which of the six
independent components of the strain are relaxed.
- True = relax to zero
- False = fixed, ignore this component
Degrees of freedom are the positions in the original undeformed
cell, plus the log of the deformation tensor (extra 3 "atoms"). This
gives forces consistent with numerical derivatives of the potential
energy with respect to the cell degrees of freedom.
For full details see:
E. B. Tadmor, G. S. Smith, N. Bernstein, and E. Kaxiras,
Phys. Rev. B 59, 235 (1999)
You can still use constraints on the atoms, e.g. FixAtoms, to
control the relaxation of the atoms.
>>> # this should be equivalent to the StrainFilter
>>> atoms = Atoms(...)
>>> atoms.set_constraint(FixAtoms(mask=[True for atom in atoms]))
>>> ecf = ExpCellFilter(atoms)
You should not attach this ExpCellFilter object to a
trajectory. Instead, create a trajectory for the atoms, and
attach it to an optimizer like this:
>>> atoms = Atoms(...)
>>> ecf = ExpCellFilter(atoms)
>>> qn = QuasiNewton(ecf)
>>> traj = Trajectory('TiO2.traj', 'w', atoms)
>>> qn.attach(traj)
>>> qn.run(fmax=0.05)
Helpful conversion table:
- 0.05 eV/A^3 = 8 GPA
- 0.003 eV/A^3 = 0.48 GPa
- 0.0006 eV/A^3 = 0.096 GPa
- 0.0003 eV/A^3 = 0.048 GPa
- 0.0001 eV/A^3 = 0.02 GPa
Additional optional arguments:
cell_factor: (DEPRECATED)
Retained for backwards compatibility, but no longer used.
hydrostatic_strain: bool (default False)
Constrain the cell by only allowing hydrostatic deformation.
The virial tensor is replaced by np.diag([np.trace(virial)]*3).
constant_volume: bool (default False)
Project out the diagonal elements of the virial tensor to allow
relaxations at constant volume, e.g. for mapping out an
energy-volume curve.
scalar_pressure: float (default 0.0)
Applied pressure to use for enthalpy pV term. As above, this
breaks energy/force consistency.
Implementation details:
The implementation is based on that of Christoph Ortner in JuLIP.jl:
https://github.com/libAtoms/JuLIP.jl/blob/expcell/src/Constraints.jl#L244
We decompose the deformation gradient as
F = exp(U) F0
x = F * F0^{-1} z = exp(U) z
If we write the energy as a function of U we can transform the
stress associated with a perturbation V into a derivative using a
linear map V -> L(U, V).
\phi( exp(U+tV) (z+tv) ) ~ \phi'(x) . (exp(U) v) + \phi'(x) .
( L(U, V) exp(-U) exp(U) z )
where
\nabla E(U) : V = [S exp(-U)'] : L(U,V)
= L'(U, S exp(-U)') : V
= L(U', S exp(-U)') : V
= L(U, S exp(-U)) : V (provided U = U')
where the : operator represents double contraction,
i.e. A:B = trace(A'B), and
F = deformation tensor - 3x3 matrix
F0 = reference deformation tensor - 3x3 matrix, np.eye(3) here
U = cell degrees of freedom used here - 3x3 matrix
V = perturbation to cell DoFs - 3x3 matrix
v = perturbation to position DoFs
x = atomic positions in deformed cell
z = atomic positions in original cell
\phi = potential energy
S = stress tensor [3x3 matrix]
L(U, V) = directional derivative of exp at U in direction V, i.e
d/dt exp(U + t V)|_{t=0} = L(U, V)
This means we can write
d/dt E(U + t V)|_{t=0} = L(U, S exp (-U)) : V
and therefore the contribution to the gradient of the energy is
\nabla E(U) / \nabla U_ij = [L(U, S exp(-U))]_ij
.. deprecated:: 3.23.0
Use :class:`~ase.filters.FrechetCellFilter` for better convergence
w.r.t. cell variables.
"""
Filter.__init__(self, atoms=atoms, indices=range(len(atoms)))
UnitCellFilter.__init__(self, atoms=atoms, mask=mask,
cell_factor=cell_factor,
hydrostatic_strain=hydrostatic_strain,
constant_volume=constant_volume,
scalar_pressure=scalar_pressure)
if cell_factor is not None:
# cell_factor used in UnitCellFilter does not affect on gradients of
# ExpCellFilter.
warn("cell_factor is deprecated")
self.cell_factor = 1.0
# We defer the scipy import to avoid high immediate import overhead
from scipy.linalg import expm, logm
self.expm = expm
self.logm = logm
def get_forces(self, **kwargs):
forces = UnitCellFilter.get_forces(self, **kwargs)
# forces on atoms are same as UnitCellFilter, we just
# need to modify the stress contribution
stress = self.atoms.get_stress(**kwargs)
volume = self.atoms.get_volume()
virial = -volume * (voigt_6_to_full_3x3_stress(stress) +
np.diag([self.scalar_pressure] * 3))
cur_deform_grad = self.deform_grad()
cur_deform_grad_log = self.logm(cur_deform_grad)
if self.hydrostatic_strain:
vtr = virial.trace()
virial = np.diag([vtr / 3.0, vtr / 3.0, vtr / 3.0])
# Zero out components corresponding to fixed lattice elements
if (self.mask != 1.0).any():
virial *= self.mask
deform_grad_log_force_naive = virial.copy()
Y = np.zeros((6, 6))
Y[0:3, 0:3] = cur_deform_grad_log
Y[3:6, 3:6] = cur_deform_grad_log
Y[0:3, 3:6] = - virial @ self.expm(-cur_deform_grad_log)
deform_grad_log_force = -self.expm(Y)[0:3, 3:6]
for (i1, i2) in [(0, 1), (0, 2), (1, 2)]:
ff = 0.5 * (deform_grad_log_force[i1, i2] +
deform_grad_log_force[i2, i1])
deform_grad_log_force[i1, i2] = ff
deform_grad_log_force[i2, i1] = ff
# check for reasonable alignment between naive and
# exact search directions
all_are_equal = np.all(np.isclose(deform_grad_log_force,
deform_grad_log_force_naive))
if all_are_equal or \
(np.sum(deform_grad_log_force * deform_grad_log_force_naive) /
np.sqrt(np.sum(deform_grad_log_force**2) *
np.sum(deform_grad_log_force_naive**2)) > 0.8):
deform_grad_log_force = deform_grad_log_force_naive
# Cauchy stress used for convergence testing
convergence_crit_stress = -(virial / volume)
if self.constant_volume:
# apply constraint to force
dglf_trace = deform_grad_log_force.trace()
np.fill_diagonal(deform_grad_log_force,
np.diag(deform_grad_log_force) - dglf_trace / 3.0)
# apply constraint to Cauchy stress used for convergence testing
ccs_trace = convergence_crit_stress.trace()
np.fill_diagonal(convergence_crit_stress,
np.diag(convergence_crit_stress) - ccs_trace / 3.0)
# pack gradients into vector
natoms = len(self.atoms)
forces[natoms:] = deform_grad_log_force
self.stress = full_3x3_to_voigt_6_stress(convergence_crit_stress)
return forces