Coverage for /builds/ase/ase/ase/calculators/tersoff.py: 95.28%
233 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1from dataclasses import dataclass
2from pathlib import Path
3from typing import Dict, List, Tuple, Type, Union
5import numpy as np
7from ase.calculators.calculator import Calculator, all_changes
8from ase.neighborlist import NeighborList
9from ase.stress import full_3x3_to_voigt_6_stress
11__author__ = 'Stefan Bringuier <stefanbringuier@gmail.com>'
12__description__ = 'LAMMPS-style native Tersoff potential for ASE'
14# Maximum/minimum exponents for numerical stability
15# in bond order calculation
16_MAX_EXP_ARG = 69.0776e0
17_MIN_EXP_ARG = -69.0776e0
20@dataclass
21class TersoffParameters:
22 """Parameters for 3 element Tersoff potential interaction.
24 Can be instantiated with either positional or keyword arguments:
25 TersoffParameters(1.0, 2.0, ...) or
26 TersoffParameters(m=1.0, gamma=2.0, ...)
27 """
29 m: float
30 gamma: float
31 lambda3: float
32 c: float
33 d: float
34 h: float
35 n: float
36 beta: float
37 lambda2: float
38 B: float
39 R: float
40 D: float
41 lambda1: float
42 A: float
44 @classmethod
45 def from_list(cls, params: List[float]) -> 'TersoffParameters':
46 """Create TersoffParameters from a list of 14 parameter values."""
47 if len(params) != 14:
48 raise ValueError(f'Expected 14 parameters, got {len(params)}')
49 return cls(*map(float, params))
52class Tersoff(Calculator):
53 """ASE Calculator for Tersoff interatomic potential.
55 .. versionadded:: 3.25.0
56 """
58 implemented_properties = [
59 'free_energy',
60 'energy',
61 'energies',
62 'forces',
63 'stress',
64 ]
66 def __init__(
67 self,
68 parameters: Dict[Tuple[str, str, str], TersoffParameters],
69 skin: float = 0.3,
70 **kwargs,
71 ) -> None:
72 """
73 Parameters
74 ----------
75 parameters : dict
76 Mapping element combinations to TersoffParameters objects::
78 {
79 ('A', 'B', 'C'): TersoffParameters(
80 m, gamma, lambda3, c, d, h, n,
81 beta, lambda2, B, R, D, lambda1, A),
82 ...
83 }
85 where ('A', 'B', 'C') are the elements involved in the interaction.
86 skin : float, default 0.3
87 The skin distance for neighbor list calculations.
88 **kwargs : dict
89 Additional parameters to be passed to
90 :class:`~ase.calculators.Calculator`.
92 """
93 Calculator.__init__(self, **kwargs)
94 self.cutoff_skin = skin
95 self.parameters = parameters
97 @classmethod
98 def from_lammps(
99 cls: Type['Tersoff'],
100 potential_file: Union[str, Path],
101 skin: float = 0.3,
102 **kwargs,
103 ) -> 'Tersoff':
104 """Make :class:`Tersoff` from a LAMMPS-style Tersoff potential file.
106 Parameters
107 ----------
108 potential_file : str or Path
109 The path to a LAMMPS-style Tersoff potential file.
110 skin : float, default 0.3
111 The skin distance for neighbor list calculations.
112 **kwargs : dict
113 Additional parameters to be passed to the
114 ASE Calculator constructor.
116 Returns
117 -------
118 :class:`Tersoff`
119 Initialized Tersoff calculator with parameters from the file.
121 """
122 parameters = cls.read_lammps_format(potential_file)
123 return cls(parameters=parameters, skin=skin, **kwargs)
125 @staticmethod
126 def read_lammps_format(
127 potential_file: Union[str, Path],
128 ) -> Dict[Tuple[str, str, str], TersoffParameters]:
129 """Read the Tersoff potential parameters from a LAMMPS-style file.
131 Parameters
132 ----------
133 potential_file : str or Path
134 Path to the LAMMPS-style Tersoff potential file
136 Returns
137 -------
138 dict
139 Dictionary mapping element combinations to TersoffParameters objects
141 """
142 block_size = 17
143 with Path(potential_file).open('r', encoding='utf-8') as fd:
144 content = (
145 ''.join(
146 [line for line in fd if not line.strip().startswith('#')]
147 )
148 .replace('\n', ' ')
149 .split()
150 )
152 if len(content) % block_size != 0:
153 raise ValueError(
154 'The potential file does not have the correct LAMMPS format.'
155 )
157 parameters: Dict[Tuple[str, str, str], TersoffParameters] = {}
158 for i in range(0, len(content), block_size):
159 block = content[i : i + block_size]
160 e1, e2, e3 = block[0], block[1], block[2]
161 current_elements = (e1, e2, e3)
162 params = map(float, block[3:])
163 parameters[current_elements] = TersoffParameters(*params)
165 return parameters
167 def set_parameters(
168 self,
169 key: Tuple[str, str, str],
170 params: TersoffParameters = None,
171 **kwargs,
172 ) -> None:
173 """Update parameters for a specific element combination.
175 Parameters
176 ----------
177 key: Tuple[str, str, str]
178 The element combination key of the parameters to be updated
179 params: TersoffParameters, optional
180 A TersoffParameters instance to completely replace the parameters
181 **kwargs:
182 Individual parameter values to update, e.g. R=2.9
184 """
185 if key not in self.parameters:
186 raise KeyError(f"Key '{key}' not found in parameters.")
188 if params is not None:
189 if kwargs:
190 raise ValueError('Cannot provide both params and kwargs.')
191 self.parameters[key] = params
192 else:
193 for name, value in kwargs.items():
194 if not hasattr(self.parameters[key], name):
195 raise ValueError(f'Invalid parameter name: {name}')
196 setattr(self.parameters[key], name, value)
198 def _update_nl(self, atoms) -> None:
199 """Update the neighbor list with the parameter R+D cutoffs.
201 Parameters
202 ----------
203 atoms: ase.Atoms
204 The atoms to calculate the neighbor list for.
206 Notes
207 -----
208 The cutoffs are determined by the parameters of the Tersoff potential.
209 Each atom's cutoff is based on the R+D values from the parameter set
210 where that atom's element appears first in the key tuple.
212 """
213 # Get cutoff for each atom based on its element type
214 cutoffs = []
216 for symbol in atoms.symbols:
217 # Find first parameter set, element is the first slot
218 param_key = next(
219 key for key in self.parameters.keys() if key[0] == symbol
220 )
221 params = self.parameters[param_key]
222 cutoffs.append(params.R + params.D)
224 self.nl = NeighborList(
225 cutoffs,
226 skin=self.cutoff_skin,
227 self_interaction=False,
228 bothways=True,
229 )
231 self.nl.update(atoms)
233 def calculate(
234 self,
235 atoms=None,
236 properties=None,
237 system_changes=all_changes,
238 ) -> None:
239 """Calculate energy, forces, and stress.
241 Notes
242 -----
243 The force and stress are calculated regardless if they are
244 requested, despite some additional overhead cost,
245 therefore they are always stored in the results dict.
247 """
248 Calculator.calculate(self, atoms, properties, system_changes)
250 # Rebuild neighbor list when any relevant system changes occur
251 checks = {'positions', 'numbers', 'cell', 'pbc'}
252 if any(change in checks for change in system_changes) or not hasattr(
253 self, 'nl'
254 ):
255 self._update_nl(atoms)
257 self.results = {}
258 energies = np.zeros(len(atoms))
259 forces = np.zeros((len(atoms), 3))
260 virial = np.zeros((3, 3))
262 # Duplicates atoms.get_distances() functionality, but uses
263 # neighbor list's pre-computed offsets for efficiency in a
264 # tight force-calculation loop rather than recompute MIC
265 for i in range(len(atoms)):
266 self._calc_atom_contribution(i, energies, forces, virial)
268 self.results['energies'] = energies
269 self.results['energy'] = self.results['free_energy'] = energies.sum()
270 self.results['forces'] = forces
271 # Virial to stress (i.e., eV/A^3)
272 if self.atoms.cell.rank == 3:
273 stress = virial / self.atoms.get_volume()
274 self.results['stress'] = full_3x3_to_voigt_6_stress(stress)
276 def _calc_atom_contribution(
277 self,
278 idx_i: int,
279 energies: np.ndarray,
280 forces: np.ndarray,
281 virial: np.ndarray,
282 ) -> None:
283 """Calculate the contributions of a single atom to the properties.
285 This function calculates the energy, force, and stress on atom i
286 by looking at i-j pair interactions and the modification made by
287 the bond order term bij with includes 3-body interaction i-j-k.
289 Parameters
290 ----------
291 idx_i: int
292 Index of atom i
293 energies: array_like
294 Site energies to be updated.
295 forces: array_like
296 Forces to be updated.
297 virial: array_like
298 Virial tensor to be updated.
300 """
301 indices, offsets = self.nl.get_neighbors(idx_i)
302 vectors = self.atoms.positions[indices]
303 vectors += offsets @ self.atoms.cell
304 vectors -= self.atoms.positions[idx_i]
305 distances = np.sqrt(np.add.reduce(vectors**2, axis=1))
307 type_i = self.atoms.symbols[idx_i]
308 for j, (idx_j, abs_rij, rij) in enumerate(
309 zip(indices, distances, vectors)
310 ):
311 type_j = self.atoms.symbols[idx_j]
312 key = (type_i, type_j, type_j)
313 params = self.parameters[key]
315 rij_hat = rij / abs_rij
317 fc = self._calc_fc(abs_rij, params.R, params.D)
318 if fc == 0.0:
319 continue
321 zeta = self._calc_zeta(type_i, j, indices, distances, vectors)
322 bij = self._calc_bij(zeta, params.beta, params.n)
323 bij_d = self._calc_bij_d(zeta, params.beta, params.n)
325 repulsive = params.A * np.exp(-params.lambda1 * abs_rij)
326 attractive = -params.B * np.exp(-params.lambda2 * abs_rij)
328 # distribute the pair energy evenly to be consistent with LAMMPS
329 energies[idx_i] += 0.25 * fc * (repulsive + bij * attractive)
330 energies[idx_j] += 0.25 * fc * (repulsive + bij * attractive)
332 dfc = self._calc_fc_d(abs_rij, params.R, params.D)
333 rep_deriv = -params.lambda1 * repulsive
334 att_deriv = -params.lambda2 * attractive
336 tmp = dfc * (repulsive + bij * attractive)
337 tmp += fc * (rep_deriv + bij * att_deriv)
339 # derivative with respect to the position of atom j
340 grad = 0.5 * tmp * rij_hat
342 forces[idx_i] += grad
343 forces[idx_j] -= grad
345 virial += np.outer(grad, rij)
347 for k, idx_k in enumerate(indices):
348 if k == j:
349 continue
351 type_k = self.atoms.symbols[idx_k]
352 key = (type_i, type_j, type_k)
353 params = self.parameters[key]
355 if distances[k] > params.R + params.D:
356 continue
358 rik = vectors[k]
360 dztdri, dztdrj, dztdrk = self._calc_zeta_d(rij, rik, params)
362 gradi = 0.5 * fc * bij_d * dztdri * attractive
363 gradj = 0.5 * fc * bij_d * dztdrj * attractive
364 gradk = 0.5 * fc * bij_d * dztdrk * attractive
366 forces[idx_i] -= gradi
367 forces[idx_j] -= gradj
368 forces[idx_k] -= gradk
370 virial += np.outer(gradj, rij)
371 virial += np.outer(gradk, rik)
373 def _calc_bij(self, zeta: float, beta: float, n: float) -> float:
374 """Calculate the bond order ``bij`` between atoms ``i`` and ``j``."""
375 tmp = beta * zeta
376 return (1.0 + tmp**n) ** (-1.0 / (2.0 * n))
378 def _calc_bij_d(self, zeta: float, beta: float, n: float) -> float:
379 """Calculate the derivative of ``bij`` with respect to ``zeta``."""
380 tmp = beta * zeta
381 return (
382 -0.5
383 * (1.0 + tmp**n) ** (-1.0 - (1.0 / (2.0 * n)))
384 * (beta * tmp ** (n - 1.0))
385 )
387 def _calc_zeta(
388 self,
389 type_i: str,
390 j: int,
391 neighbors: np.ndarray,
392 distances: np.ndarray,
393 vectors: np.ndarray,
394 ) -> float:
395 """Calculate ``zeta_ij``."""
396 idx_j = neighbors[j]
397 type_j = self.atoms.symbols[idx_j]
398 abs_rij = distances[j]
400 zeta = 0.0
402 for k, idx_k in enumerate(neighbors):
403 if k == j:
404 continue
406 type_k = self.atoms.symbols[idx_k]
407 key = (type_i, type_j, type_k)
408 params = self.parameters[key]
410 abs_rik = distances[k]
411 if abs_rik > params.R + params.D:
412 continue
414 costheta = np.dot(vectors[j], vectors[k]) / (abs_rij * abs_rik)
415 fc_ik = self._calc_fc(abs_rik, params.R, params.D)
417 g_theta = self._calc_gijk(costheta, params)
419 # Calculate the exponential for the bond order zeta term
420 # This is the term that modifies the bond order based
421 # on the distance between atoms i-j and i-k. Tresholds are
422 # used to prevent overflow/underflow.
423 arg = (params.lambda3 * (abs_rij - abs_rik)) ** params.m
424 if arg > _MAX_EXP_ARG:
425 ex_delr = 1.0e30
426 elif arg < _MIN_EXP_ARG:
427 ex_delr = 0.0
428 else:
429 ex_delr = np.exp(arg)
431 zeta += fc_ik * g_theta * ex_delr
433 return zeta
435 def _calc_gijk(self, costheta: float, params: TersoffParameters) -> float:
436 r"""Calculate the angular function ``g`` for the Tersoff potential.
438 .. math::
439 g(\theta) = \gamma \left( 1 + \frac{c^2}{d^2}
440 - \frac{c^2}{d^2 + (h - \cos \theta)^2} \right)
442 where :math:`\theta` is the angle between the bond vector
443 and the vector of atom i and its neighbors j-k.
444 """
445 c2 = params.c * params.c
446 d2 = params.d * params.d
447 hcth = params.h - costheta
448 return params.gamma * (1.0 + c2 / d2 - c2 / (d2 + hcth**2))
450 def _calc_gijk_d(self, costheta: float, params: TersoffParameters) -> float:
451 """Calculate the derivative of ``g`` with respect to ``costheta``."""
452 c2 = params.c * params.c
453 d2 = params.d * params.d
454 hcth = params.h - costheta
455 numerator = -2.0 * params.gamma * c2 * hcth
456 denominator = (d2 + hcth**2) ** 2
457 return numerator / denominator
459 def _calc_fc(self, r: np.floating, R: float, D: float) -> float:
460 """Calculate the cutoff function."""
461 if r > R + D:
462 return 0.0
463 if r < R - D:
464 return 1.0
465 return 0.5 * (1.0 - np.sin(np.pi * (r - R) / (2.0 * D)))
467 def _calc_fc_d(self, r: np.floating, R: float, D: float) -> float:
468 """Calculate cutoff function derivative with respect to ``r``."""
469 if r > R + D or r < R - D:
470 return 0.0
471 return -0.25 * np.pi / D * np.cos(np.pi * (r - R) / (2.0 * D))
473 def _calc_zeta_d(
474 self,
475 rij: np.ndarray,
476 rik: np.ndarray,
477 params: TersoffParameters,
478 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
479 """Calculate the derivatives of ``zeta``.
481 Returns
482 -------
483 dri : ndarray of shape (3,), dtype float
484 Derivative with respect to the position of atom ``i``.
485 drj : ndarray of shape (3,), dtype float
486 Derivative with respect to the position of atom ``j``.
487 drk : ndarray of shape (3,), dtype float
488 Derivative with respect to the position of atom ``k``.
490 """
491 lam3 = params.lambda3
492 m = params.m
494 abs_rij = np.linalg.norm(rij)
495 abs_rik = np.linalg.norm(rik)
497 rij_hat = rij / abs_rij
498 rik_hat = rik / abs_rik
500 fcik = self._calc_fc(abs_rik, params.R, params.D)
501 dfcik = self._calc_fc_d(abs_rik, params.R, params.D)
503 tmp = (lam3 * (abs_rij - abs_rik)) ** m
504 if tmp > _MAX_EXP_ARG:
505 ex_delr = 1.0e30
506 elif tmp < _MIN_EXP_ARG:
507 ex_delr = 0.0
508 else:
509 ex_delr = np.exp(tmp)
511 ex_delr_d = m * lam3**m * (abs_rij - abs_rik) ** (m - 1) * ex_delr
513 costheta = rij_hat @ rik_hat
514 gijk = self._calc_gijk(costheta, params)
515 gijk_d = self._calc_gijk_d(costheta, params)
517 dcosdri, dcosdrj, dcosdrk = self._calc_costheta_d(rij, rik)
519 dri = -dfcik * gijk * ex_delr * rik_hat
520 dri += fcik * gijk_d * ex_delr * dcosdri
521 dri += fcik * gijk * ex_delr_d * rik_hat
522 dri -= fcik * gijk * ex_delr_d * rij_hat
524 drj = fcik * gijk_d * ex_delr * dcosdrj
525 drj += fcik * gijk * ex_delr_d * rij_hat
527 drk = dfcik * gijk * ex_delr * rik_hat
528 drk += fcik * gijk_d * ex_delr * dcosdrk
529 drk -= fcik * gijk * ex_delr_d * rik_hat
531 return dri, drj, drk
533 def _calc_costheta_d(
534 self,
535 rij: np.ndarray,
536 rik: np.ndarray,
537 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
538 r"""Calculate the derivatives of ``costheta``.
540 If
542 .. math::
543 \cos \theta = \frac{\mathbf{u} \cdot \mathbf{v}}{u v}
545 Then
547 .. math::
548 \frac{\partial \cos \theta}{\partial \mathbf{u}}
549 = \frac{\mathbf{v}}{u v}
550 - \frac{\mathbf{u} \cdot \mathbf{v}}{v} \cdot \frac{\mathbf{u}}{u^3}
551 = \frac{\mathbf{v}}{u v} - \frac{\cos \theta}{u^2} \mathbf{u}
553 Parameters
554 ----------
555 rij : ndarray of shape (3,), dtype float
556 Vector from atoms ``i`` to ``j``.
557 rik : ndarray of shape (3,), dtype float
558 Vector from atoms ``i`` to ``k``.
560 Returns
561 -------
562 dri : ndarray of shape (3,), dtype float
563 Derivative with respect to the position of atom ``i``.
564 drj : ndarray of shape (3,), dtype float
565 Derivative with respect to the position of atom ``j``.
566 drk : ndarray of shape (3,), dtype float
567 Derivative with respect to the position of atom ``k``.
569 """
570 abs_rij = np.linalg.norm(rij)
571 abs_rik = np.linalg.norm(rik)
572 costheta = (rij @ rik) / (abs_rij * abs_rik)
573 drj = (rik / abs_rik - costheta * rij / abs_rij) / abs_rij
574 drk = (rij / abs_rij - costheta * rik / abs_rik) / abs_rik
575 dri = -(drj + drk)
576 return dri, drj, drk