Coverage for ase / calculators / tersoff.py: 95.26%
232 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1from dataclasses import dataclass
2from pathlib import Path
4import numpy as np
6from ase.calculators.calculator import Calculator, all_changes
7from ase.neighborlist import NeighborList
8from ase.stress import full_3x3_to_voigt_6_stress
10__author__ = 'Stefan Bringuier <stefanbringuier@gmail.com>'
11__description__ = 'LAMMPS-style native Tersoff potential for ASE'
13# Maximum/minimum exponents for numerical stability
14# in bond order calculation
15_MAX_EXP_ARG = 69.0776e0
16_MIN_EXP_ARG = -69.0776e0
19@dataclass
20class TersoffParameters:
21 """Parameters for 3 element Tersoff potential interaction.
23 Can be instantiated with either positional or keyword arguments:
24 TersoffParameters(1.0, 2.0, ...) or
25 TersoffParameters(m=1.0, gamma=2.0, ...)
26 """
28 m: float
29 gamma: float
30 lambda3: float
31 c: float
32 d: float
33 h: float
34 n: float
35 beta: float
36 lambda2: float
37 B: float
38 R: float
39 D: float
40 lambda1: float
41 A: float
43 @classmethod
44 def from_list(cls, params: list[float]) -> 'TersoffParameters':
45 """Create TersoffParameters from a list of 14 parameter values."""
46 if len(params) != 14:
47 raise ValueError(f'Expected 14 parameters, got {len(params)}')
48 return cls(*map(float, params))
51class Tersoff(Calculator):
52 """ASE Calculator for Tersoff interatomic potential.
54 .. versionadded:: 3.25.0
55 """
57 implemented_properties = [
58 'free_energy',
59 'energy',
60 'energies',
61 'forces',
62 'stress',
63 ]
65 def __init__(
66 self,
67 parameters: dict[tuple[str, str, str], TersoffParameters],
68 skin: float = 0.3,
69 **kwargs,
70 ) -> None:
71 """
72 Parameters
73 ----------
74 parameters : dict
75 Mapping element combinations to TersoffParameters objects::
77 {
78 ('A', 'B', 'C'): TersoffParameters(
79 m, gamma, lambda3, c, d, h, n,
80 beta, lambda2, B, R, D, lambda1, A),
81 ...
82 }
84 where ('A', 'B', 'C') are the elements involved in the interaction.
85 skin : float, default 0.3
86 The skin distance for neighbor list calculations.
87 **kwargs : dict
88 Additional parameters to be passed to
89 :class:`~ase.calculators.Calculator`.
91 """
92 Calculator.__init__(self, **kwargs)
93 self.cutoff_skin = skin
94 self.parameters = parameters
96 @classmethod
97 def from_lammps(
98 cls: type['Tersoff'],
99 potential_file: str | Path,
100 skin: float = 0.3,
101 **kwargs,
102 ) -> 'Tersoff':
103 """Make :class:`Tersoff` from a LAMMPS-style Tersoff potential file.
105 Parameters
106 ----------
107 potential_file : str or Path
108 The path to a LAMMPS-style Tersoff potential file.
109 skin : float, default 0.3
110 The skin distance for neighbor list calculations.
111 **kwargs : dict
112 Additional parameters to be passed to the
113 ASE Calculator constructor.
115 Returns
116 -------
117 :class:`Tersoff`
118 Initialized Tersoff calculator with parameters from the file.
120 """
121 parameters = cls.read_lammps_format(potential_file)
122 return cls(parameters=parameters, skin=skin, **kwargs)
124 @staticmethod
125 def read_lammps_format(
126 potential_file: str | Path,
127 ) -> dict[tuple[str, str, str], TersoffParameters]:
128 """Read the Tersoff potential parameters from a LAMMPS-style file.
130 Parameters
131 ----------
132 potential_file : str or Path
133 Path to the LAMMPS-style Tersoff potential file
135 Returns
136 -------
137 dict
138 Dictionary mapping element combinations to TersoffParameters objects
140 """
141 block_size = 17
142 with Path(potential_file).open('r', encoding='utf-8') as fd:
143 content = (
144 ''.join(
145 [line for line in fd if not line.strip().startswith('#')]
146 )
147 .replace('\n', ' ')
148 .split()
149 )
151 if len(content) % block_size != 0:
152 raise ValueError(
153 'The potential file does not have the correct LAMMPS format.'
154 )
156 parameters: dict[tuple[str, str, str], TersoffParameters] = {}
157 for i in range(0, len(content), block_size):
158 block = content[i : i + block_size]
159 e1, e2, e3 = block[0], block[1], block[2]
160 current_elements = (e1, e2, e3)
161 params = map(float, block[3:])
162 parameters[current_elements] = TersoffParameters(*params)
164 return parameters
166 def set_parameters(
167 self,
168 key: tuple[str, str, str],
169 params: TersoffParameters | None = None,
170 **kwargs,
171 ) -> None:
172 """Update parameters for a specific element combination.
174 Parameters
175 ----------
176 key: Tuple[str, str, str]
177 The element combination key of the parameters to be updated
178 params: TersoffParameters, optional
179 A TersoffParameters instance to completely replace the parameters
180 **kwargs:
181 Individual parameter values to update, e.g. R=2.9
183 """
184 if key not in self.parameters:
185 raise KeyError(f"Key '{key}' not found in parameters.")
187 if params is not None:
188 if kwargs:
189 raise ValueError('Cannot provide both params and kwargs.')
190 self.parameters[key] = params
191 else:
192 for name, value in kwargs.items():
193 if not hasattr(self.parameters[key], name):
194 raise ValueError(f'Invalid parameter name: {name}')
195 setattr(self.parameters[key], name, value)
197 def _update_nl(self, atoms) -> None:
198 """Update the neighbor list with the parameter R+D cutoffs.
200 Parameters
201 ----------
202 atoms: ase.Atoms
203 The atoms to calculate the neighbor list for.
205 Notes
206 -----
207 The cutoffs are determined by the parameters of the Tersoff potential.
208 Each atom's cutoff is based on the R+D values from the parameter set
209 where that atom's element appears first in the key tuple.
211 """
212 # Get cutoff for each atom based on its element type
213 cutoffs = []
215 for symbol in atoms.symbols:
216 # Find first parameter set, element is the first slot
217 param_key = next(
218 key for key in self.parameters.keys() if key[0] == symbol
219 )
220 params = self.parameters[param_key]
221 cutoffs.append(params.R + params.D)
223 self.nl = NeighborList(
224 cutoffs,
225 skin=self.cutoff_skin,
226 self_interaction=False,
227 bothways=True,
228 )
230 self.nl.update(atoms)
232 def calculate(
233 self,
234 atoms=None,
235 properties=None,
236 system_changes=all_changes,
237 ) -> None:
238 """Calculate energy, forces, and stress.
240 Notes
241 -----
242 The force and stress are calculated regardless if they are
243 requested, despite some additional overhead cost,
244 therefore they are always stored in the results dict.
246 """
247 Calculator.calculate(self, atoms, properties, system_changes)
249 # Rebuild neighbor list when any relevant system changes occur
250 checks = {'positions', 'numbers', 'cell', 'pbc'}
251 if any(change in checks for change in system_changes) or not hasattr(
252 self, 'nl'
253 ):
254 self._update_nl(atoms)
256 self.results = {}
257 energies = np.zeros(len(atoms))
258 forces = np.zeros((len(atoms), 3))
259 virial = np.zeros((3, 3))
261 # Duplicates atoms.get_distances() functionality, but uses
262 # neighbor list's pre-computed offsets for efficiency in a
263 # tight force-calculation loop rather than recompute MIC
264 for i in range(len(atoms)):
265 self._calc_atom_contribution(i, energies, forces, virial)
267 self.results['energies'] = energies
268 self.results['energy'] = self.results['free_energy'] = energies.sum()
269 self.results['forces'] = forces
270 # Virial to stress (i.e., eV/A^3)
271 if self.atoms.cell.rank == 3:
272 stress = virial / self.atoms.get_volume()
273 self.results['stress'] = full_3x3_to_voigt_6_stress(stress)
275 def _calc_atom_contribution(
276 self,
277 idx_i: int,
278 energies: np.ndarray,
279 forces: np.ndarray,
280 virial: np.ndarray,
281 ) -> None:
282 """Calculate the contributions of a single atom to the properties.
284 This function calculates the energy, force, and stress on atom i
285 by looking at i-j pair interactions and the modification made by
286 the bond order term bij with includes 3-body interaction i-j-k.
288 Parameters
289 ----------
290 idx_i: int
291 Index of atom i
292 energies: array_like
293 Site energies to be updated.
294 forces: array_like
295 Forces to be updated.
296 virial: array_like
297 Virial tensor to be updated.
299 """
300 indices, offsets = self.nl.get_neighbors(idx_i)
301 vectors = self.atoms.positions[indices]
302 vectors += offsets @ self.atoms.cell
303 vectors -= self.atoms.positions[idx_i]
304 distances = np.sqrt(np.add.reduce(vectors**2, axis=1))
306 type_i = self.atoms.symbols[idx_i]
307 for j, (idx_j, abs_rij, rij) in enumerate(
308 zip(indices, distances, vectors)
309 ):
310 type_j = self.atoms.symbols[idx_j]
311 key = (type_i, type_j, type_j)
312 params = self.parameters[key]
314 rij_hat = rij / abs_rij
316 fc = self._calc_fc(abs_rij, params.R, params.D)
317 if fc == 0.0:
318 continue
320 zeta = self._calc_zeta(type_i, j, indices, distances, vectors)
321 bij = self._calc_bij(zeta, params.beta, params.n)
322 bij_d = self._calc_bij_d(zeta, params.beta, params.n)
324 repulsive = params.A * np.exp(-params.lambda1 * abs_rij)
325 attractive = -params.B * np.exp(-params.lambda2 * abs_rij)
327 # distribute the pair energy evenly to be consistent with LAMMPS
328 energies[idx_i] += 0.25 * fc * (repulsive + bij * attractive)
329 energies[idx_j] += 0.25 * fc * (repulsive + bij * attractive)
331 dfc = self._calc_fc_d(abs_rij, params.R, params.D)
332 rep_deriv = -params.lambda1 * repulsive
333 att_deriv = -params.lambda2 * attractive
335 tmp = dfc * (repulsive + bij * attractive)
336 tmp += fc * (rep_deriv + bij * att_deriv)
338 # derivative with respect to the position of atom j
339 grad = 0.5 * tmp * rij_hat
341 forces[idx_i] += grad
342 forces[idx_j] -= grad
344 virial += np.outer(grad, rij)
346 for k, idx_k in enumerate(indices):
347 if k == j:
348 continue
350 type_k = self.atoms.symbols[idx_k]
351 key = (type_i, type_j, type_k)
352 params = self.parameters[key]
354 if distances[k] > params.R + params.D:
355 continue
357 rik = vectors[k]
359 dztdri, dztdrj, dztdrk = self._calc_zeta_d(rij, rik, params)
361 gradi = 0.5 * fc * bij_d * dztdri * attractive
362 gradj = 0.5 * fc * bij_d * dztdrj * attractive
363 gradk = 0.5 * fc * bij_d * dztdrk * attractive
365 forces[idx_i] -= gradi
366 forces[idx_j] -= gradj
367 forces[idx_k] -= gradk
369 virial += np.outer(gradj, rij)
370 virial += np.outer(gradk, rik)
372 def _calc_bij(self, zeta: float, beta: float, n: float) -> float:
373 """Calculate the bond order ``bij`` between atoms ``i`` and ``j``."""
374 tmp = beta * zeta
375 return (1.0 + tmp**n) ** (-1.0 / (2.0 * n))
377 def _calc_bij_d(self, zeta: float, beta: float, n: float) -> float:
378 """Calculate the derivative of ``bij`` with respect to ``zeta``."""
379 tmp = beta * zeta
380 return (
381 -0.5
382 * (1.0 + tmp**n) ** (-1.0 - (1.0 / (2.0 * n)))
383 * (beta * tmp ** (n - 1.0))
384 )
386 def _calc_zeta(
387 self,
388 type_i: str,
389 j: int,
390 neighbors: np.ndarray,
391 distances: np.ndarray,
392 vectors: np.ndarray,
393 ) -> float:
394 """Calculate ``zeta_ij``."""
395 idx_j = neighbors[j]
396 type_j = self.atoms.symbols[idx_j]
397 abs_rij = distances[j]
399 zeta = 0.0
401 for k, idx_k in enumerate(neighbors):
402 if k == j:
403 continue
405 type_k = self.atoms.symbols[idx_k]
406 key = (type_i, type_j, type_k)
407 params = self.parameters[key]
409 abs_rik = distances[k]
410 if abs_rik > params.R + params.D:
411 continue
413 costheta = np.dot(vectors[j], vectors[k]) / (abs_rij * abs_rik)
414 fc_ik = self._calc_fc(abs_rik, params.R, params.D)
416 g_theta = self._calc_gijk(costheta, params)
418 # Calculate the exponential for the bond order zeta term
419 # This is the term that modifies the bond order based
420 # on the distance between atoms i-j and i-k. Tresholds are
421 # used to prevent overflow/underflow.
422 arg = (params.lambda3 * (abs_rij - abs_rik)) ** params.m
423 if arg > _MAX_EXP_ARG:
424 ex_delr = 1.0e30
425 elif arg < _MIN_EXP_ARG:
426 ex_delr = 0.0
427 else:
428 ex_delr = np.exp(arg)
430 zeta += fc_ik * g_theta * ex_delr
432 return zeta
434 def _calc_gijk(self, costheta: float, params: TersoffParameters) -> float:
435 r"""Calculate the angular function ``g`` for the Tersoff potential.
437 .. math::
438 g(\theta) = \gamma \left( 1 + \frac{c^2}{d^2}
439 - \frac{c^2}{d^2 + (h - \cos \theta)^2} \right)
441 where :math:`\theta` is the angle between the bond vector
442 and the vector of atom i and its neighbors j-k.
443 """
444 c2 = params.c * params.c
445 d2 = params.d * params.d
446 hcth = params.h - costheta
447 return params.gamma * (1.0 + c2 / d2 - c2 / (d2 + hcth**2))
449 def _calc_gijk_d(self, costheta: float, params: TersoffParameters) -> float:
450 """Calculate the derivative of ``g`` with respect to ``costheta``."""
451 c2 = params.c * params.c
452 d2 = params.d * params.d
453 hcth = params.h - costheta
454 numerator = -2.0 * params.gamma * c2 * hcth
455 denominator = (d2 + hcth**2) ** 2
456 return numerator / denominator
458 def _calc_fc(self, r: np.floating, R: float, D: float) -> float:
459 """Calculate the cutoff function."""
460 if r > R + D:
461 return 0.0
462 if r < R - D:
463 return 1.0
464 return 0.5 * (1.0 - np.sin(np.pi * (r - R) / (2.0 * D)))
466 def _calc_fc_d(self, r: np.floating, R: float, D: float) -> float:
467 """Calculate cutoff function derivative with respect to ``r``."""
468 if r > R + D or r < R - D:
469 return 0.0
470 return -0.25 * np.pi / D * np.cos(np.pi * (r - R) / (2.0 * D))
472 def _calc_zeta_d(
473 self,
474 rij: np.ndarray,
475 rik: np.ndarray,
476 params: TersoffParameters,
477 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
478 """Calculate the derivatives of ``zeta``.
480 Returns
481 -------
482 dri : ndarray of shape (3,), dtype float
483 Derivative with respect to the position of atom ``i``.
484 drj : ndarray of shape (3,), dtype float
485 Derivative with respect to the position of atom ``j``.
486 drk : ndarray of shape (3,), dtype float
487 Derivative with respect to the position of atom ``k``.
489 """
490 lam3 = params.lambda3
491 m = params.m
493 abs_rij = np.linalg.norm(rij)
494 abs_rik = np.linalg.norm(rik)
496 rij_hat = rij / abs_rij
497 rik_hat = rik / abs_rik
499 fcik = self._calc_fc(abs_rik, params.R, params.D)
500 dfcik = self._calc_fc_d(abs_rik, params.R, params.D)
502 tmp = (lam3 * (abs_rij - abs_rik)) ** m
503 if tmp > _MAX_EXP_ARG:
504 ex_delr = 1.0e30
505 elif tmp < _MIN_EXP_ARG:
506 ex_delr = 0.0
507 else:
508 ex_delr = np.exp(tmp)
510 ex_delr_d = m * lam3**m * (abs_rij - abs_rik) ** (m - 1) * ex_delr
512 costheta = rij_hat @ rik_hat
513 gijk = self._calc_gijk(costheta, params)
514 gijk_d = self._calc_gijk_d(costheta, params)
516 dcosdri, dcosdrj, dcosdrk = self._calc_costheta_d(rij, rik)
518 dri = -dfcik * gijk * ex_delr * rik_hat
519 dri += fcik * gijk_d * ex_delr * dcosdri
520 dri += fcik * gijk * ex_delr_d * rik_hat
521 dri -= fcik * gijk * ex_delr_d * rij_hat
523 drj = fcik * gijk_d * ex_delr * dcosdrj
524 drj += fcik * gijk * ex_delr_d * rij_hat
526 drk = dfcik * gijk * ex_delr * rik_hat
527 drk += fcik * gijk_d * ex_delr * dcosdrk
528 drk -= fcik * gijk * ex_delr_d * rik_hat
530 return dri, drj, drk
532 def _calc_costheta_d(
533 self,
534 rij: np.ndarray,
535 rik: np.ndarray,
536 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
537 r"""Calculate the derivatives of ``costheta``.
539 If
541 .. math::
542 \cos \theta = \frac{\mathbf{u} \cdot \mathbf{v}}{u v}
544 Then
546 .. math::
547 \frac{\partial \cos \theta}{\partial \mathbf{u}}
548 = \frac{\mathbf{v}}{u v}
549 - \frac{\mathbf{u} \cdot \mathbf{v}}{v} \cdot \frac{\mathbf{u}}{u^3}
550 = \frac{\mathbf{v}}{u v} - \frac{\cos \theta}{u^2} \mathbf{u}
552 Parameters
553 ----------
554 rij : ndarray of shape (3,), dtype float
555 Vector from atoms ``i`` to ``j``.
556 rik : ndarray of shape (3,), dtype float
557 Vector from atoms ``i`` to ``k``.
559 Returns
560 -------
561 dri : ndarray of shape (3,), dtype float
562 Derivative with respect to the position of atom ``i``.
563 drj : ndarray of shape (3,), dtype float
564 Derivative with respect to the position of atom ``j``.
565 drk : ndarray of shape (3,), dtype float
566 Derivative with respect to the position of atom ``k``.
568 """
569 abs_rij = np.linalg.norm(rij)
570 abs_rik = np.linalg.norm(rik)
571 costheta = (rij @ rik) / (abs_rij * abs_rik)
572 drj = (rik / abs_rik - costheta * rij / abs_rij) / abs_rij
573 drk = (rij / abs_rij - costheta * rik / abs_rik) / abs_rik
574 dri = -(drj + drk)
575 return dri, drj, drk