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