Coverage for ase / calculators / fd.py: 100.00%
70 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1"""Module for `FiniteDifferenceCalculator`."""
3from collections.abc import Iterable
4from functools import partial
6import numpy as np
8from ase import Atoms
9from ase.calculators.calculator import BaseCalculator, all_properties
12class FiniteDifferenceCalculator(BaseCalculator):
13 """Wrapper calculator using the finite-difference method.
15 The forces and the stress are computed using the finite-difference method.
17 .. versionadded:: 3.24.0
18 """
20 implemented_properties = all_properties
22 def __init__(
23 self,
24 calc: BaseCalculator,
25 eps_disp: float | None = 1e-6,
26 eps_strain: float | None = 1e-6,
27 *,
28 force_consistent: bool = True,
29 ) -> None:
30 """
32 Parameters
33 ----------
34 calc : :class:`~ase.calculators.calculator.BaseCalculator`
35 ASE Calculator object to be wrapped.
36 eps_disp : Optional[float], default 1e-6
37 Displacement used for computing forces.
38 If :py:obj:`None`, analytical forces are computed.
39 eps_strain : Optional[float], default 1e-6
40 Strain used for computing stress.
41 If :py:obj:`None`, analytical stress is computed.
42 force_consistent : bool, default :py:obj:`True`
43 If :py:obj:`True`, the energies consistent with the forces are used
44 for finite-difference calculations.
46 """
47 super().__init__()
48 self.calc = calc
49 self.eps_disp = eps_disp
50 self.eps_strain = eps_strain
51 self.force_consistent = force_consistent
53 def calculate(self, atoms: Atoms, properties, system_changes) -> None:
54 atoms = atoms.copy() # copy to not mess up original `atoms`
55 atoms.calc = self.calc
56 self.results = {}
57 self.results['energy'] = self.calc.get_potential_energy(atoms)
58 for key in ['free_energy']:
59 if key in self.calc.results:
60 self.results[key] = self.calc.results[key]
61 if self.eps_disp is None:
62 self.results['forces'] = self.calc.get_forces(atoms)
63 else:
64 self.results['forces'] = calculate_numerical_forces(
65 atoms,
66 eps=self.eps_disp,
67 force_consistent=self.force_consistent,
68 )
69 if self.eps_strain is None:
70 self.results['stress'] = self.calc.get_stress(atoms)
71 else:
72 self.results['stress'] = calculate_numerical_stress(
73 atoms,
74 eps=self.eps_strain,
75 force_consistent=self.force_consistent,
76 )
79def _numeric_force(
80 atoms: Atoms,
81 iatom: int,
82 icart: int,
83 eps: float = 1e-6,
84 *,
85 force_consistent: bool = False,
86) -> float:
87 """Calculate numerical force on a specific atom along a specific direction.
89 Parameters
90 ----------
91 atoms : :class:`~ase.Atoms`
92 ASE :class:`~ase.Atoms` object.
93 iatom : int
94 Index of atoms.
95 icart : {0, 1, 2}
96 Index of Cartesian component.
97 eps : float, default 1e-6
98 Displacement.
99 force_consistent : bool, default :py:obj:`False`
100 If :py:obj:`True`, the energies consistent with the forces are used for
101 finite-difference calculations.
103 """
104 p0 = atoms.get_positions()
105 p = p0.copy()
106 p[iatom, icart] = p0[iatom, icart] + eps
107 atoms.set_positions(p, apply_constraint=False)
108 eplus = atoms.get_potential_energy(force_consistent=force_consistent)
109 p[iatom, icart] = p0[iatom, icart] - eps
110 atoms.set_positions(p, apply_constraint=False)
111 eminus = atoms.get_potential_energy(force_consistent=force_consistent)
112 atoms.set_positions(p0, apply_constraint=False)
113 return (eminus - eplus) / (2 * eps)
116def calculate_numerical_forces(
117 atoms: Atoms,
118 eps: float = 1e-6,
119 iatoms: Iterable[int] | None = None,
120 icarts: Iterable[int] | None = None,
121 *,
122 force_consistent: bool = False,
123) -> np.ndarray:
124 """Calculate forces numerically based on the finite-difference method.
126 Parameters
127 ----------
128 atoms : :class:`~ase.Atoms`
129 ASE :class:`~ase.Atoms` object.
130 eps : float, default 1e-6
131 Displacement.
132 iatoms : Optional[Iterable[int]]
133 Indices of atoms for which forces are computed.
134 By default, all atoms are considered.
135 icarts : Optional[Iterable[int]]
136 Indices of Cartesian coordinates for which forces are computed.
137 By default, all three coordinates are considered.
138 force_consistent : bool, default :py:obj:`False`
139 If :py:obj:`True`, the energies consistent with the forces are used for
140 finite-difference calculations.
142 Returns
143 -------
144 forces : np.ndarray
145 Forces computed numerically based on the finite-difference method.
147 """
148 if iatoms is None:
149 iatoms = range(len(atoms))
150 if icarts is None:
151 icarts = [0, 1, 2]
152 f = partial(_numeric_force, eps=eps, force_consistent=force_consistent)
153 forces = [[f(atoms, iatom, icart) for icart in icarts] for iatom in iatoms]
154 return np.array(forces)
157def calculate_numerical_stress(
158 atoms: Atoms,
159 eps: float = 1e-6,
160 voigt: bool = True,
161 *,
162 force_consistent: bool = True,
163) -> np.ndarray:
164 """Calculate stress numerically based on the finite-difference method.
166 Parameters
167 ----------
168 atoms : :class:`~ase.Atoms`
169 ASE :class:`~ase.Atoms` object.
170 eps : float, default 1e-6
171 Strain in the Voigt notation.
172 voigt : bool, default :py:obj:`True`
173 If :py:obj:`True`, the stress is returned in the Voigt notation.
174 force_consistent : bool, default :py:obj:`True`
175 If :py:obj:`True`, the energies consistent with the forces are used for
176 finite-difference calculations.
178 Returns
179 -------
180 stress : np.ndarray
181 Stress computed numerically based on the finite-difference method.
183 """
184 stress = np.zeros((3, 3), dtype=float)
186 cell = atoms.cell.copy()
187 volume = atoms.get_volume()
188 for i in range(3):
189 x = np.eye(3)
190 x[i, i] = 1.0 + eps
191 atoms.set_cell(cell @ x, scale_atoms=True)
192 eplus = atoms.get_potential_energy(force_consistent=force_consistent)
194 x[i, i] = 1.0 - eps
195 atoms.set_cell(cell @ x, scale_atoms=True)
196 eminus = atoms.get_potential_energy(force_consistent=force_consistent)
198 stress[i, i] = (eplus - eminus) / (2 * eps * volume)
199 x[i, i] = 1.0
201 j = i - 2
202 x[i, j] = x[j, i] = +0.5 * eps
203 atoms.set_cell(cell @ x, scale_atoms=True)
204 eplus = atoms.get_potential_energy(force_consistent=force_consistent)
206 x[i, j] = x[j, i] = -0.5 * eps
207 atoms.set_cell(cell @ x, scale_atoms=True)
208 eminus = atoms.get_potential_energy(force_consistent=force_consistent)
210 stress[i, j] = stress[j, i] = (eplus - eminus) / (2 * eps * volume)
212 atoms.set_cell(cell, scale_atoms=True)
214 return stress.flat[[0, 4, 8, 5, 2, 1]] if voigt else stress