Filters¶

The Filter class¶

Constraints can also be applied via filters, which acts as a wrapper around an atoms object. A typical use case will look like this:

 -------       --------       ----------
|       |     |        |     |          |
| Atoms |<----| Filter |<----| Dynamics |
|       |     |        |     |          |
 -------       --------       ----------

and in Python this would be:

>>> atoms = Atoms(...)
>>> filter = Filter(atoms, ...)
>>> dyn = Dynamics(filter, ...)

This class hides some of the atoms in an Atoms object.

class ase.filters.Filter(atoms, indices=None, mask=None)[source]¶

You must supply either the indices of the atoms that should be kept visible or a mask. The mask is a list of booleans, one for each atom, being true if the atom should be kept visible.

Example of use:

>>> from ase import Atoms, Filter
>>> atoms=Atoms(positions=[[ 0    , 0    , 0],
...                        [ 0.773, 0.600, 0],
...                        [-0.773, 0.600, 0]],
...             symbols='OH2')
>>> f1 = Filter(atoms, indices=[1, 2])
>>> f2 = Filter(atoms, mask=[0, 1, 1])
>>> f3 = Filter(atoms, mask=[a.Z == 1 for a in atoms])
>>> f1.get_positions()
[[ 0.773  0.6    0.   ]
 [-0.773  0.6    0.   ]]

In all three filters only the hydrogen atoms are made visible. When asking for the positions only the positions of the hydrogen atoms are returned.

The UnitCellFilter class¶

The unit cell filter is for optimizing positions and unit cell simultaneously. Note that ExpCellFilter will probably perform better.

class ase.filters.UnitCellFilter(atoms, mask=None, cell_factor=None, hydrostatic_strain=False, constant_volume=False, orig_cell=None, scalar_pressure=0.0)[source]¶

Modify the supercell and the atom positions.

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.

The StrainFilter class¶

The strain filter is for optimizing the unit cell while keeping scaled positions fixed.

class ase.filters.StrainFilter(atoms, mask=None, include_ideal_gas=False)[source]¶

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.

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].

The ExpCellFilter class¶

The exponential cell filter is an improved UnitCellFilter which is parameter free.

class ase.filters.ExpCellFilter(atoms, mask=None, cell_factor=None, hydrostatic_strain=False, constant_volume=False, scalar_pressure=0.0)[source]¶

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 since version 3.23.0: Use FrechetCellFilter for better convergence w.r.t. cell variables.

The FrechetCellFilter class¶

class ase.filters.FrechetCellFilter(atoms, mask=None, exp_cell_factor=None, hydrostatic_strain=False, constant_volume=False, scalar_pressure=0.0)[source]¶

Modify the supercell and the atom positions.

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