Coverage for /builds/ase/ase/ase/calculators/lj.py: 100.00%
69 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import numpy as np
5from ase.calculators.calculator import Calculator, all_changes
6from ase.neighborlist import NeighborList
7from ase.stress import full_3x3_to_voigt_6_stress
10class LennardJones(Calculator):
11 """Lennard Jones potential calculator
13 see https://en.wikipedia.org/wiki/Lennard-Jones_potential
15 The fundamental definition of this potential is a pairwise energy:
17 ``u_ij = 4 epsilon ( sigma^12/r_ij^12 - sigma^6/r_ij^6 )``
19 For convenience, we'll use d_ij to refer to "distance vector" and
20 ``r_ij`` to refer to "scalar distance". So, with position vectors `r_i`:
22 ``r_ij = | r_j - r_i | = | d_ij |``
24 Therefore:
26 ``d r_ij / d d_ij = + d_ij / r_ij``
27 ``d r_ij / d d_i = - d_ij / r_ij``
29 The derivative of u_ij is:
31 ::
33 d u_ij / d r_ij
34 = (-24 epsilon / r_ij) ( 2 sigma^12/r_ij^12 - sigma^6/r_ij^6 )
36 We can define a "pairwise force"
38 ``f_ij = d u_ij / d d_ij = d u_ij / d r_ij * d_ij / r_ij``
40 The terms in front of d_ij are combined into a "general derivative".
42 ``du_ij = (d u_ij / d d_ij) / r_ij``
44 We do this for convenience: `du_ij` is purely scalar The pairwise force is:
46 ``f_ij = du_ij * d_ij``
48 The total force on an atom is:
50 ``f_i = sum_(j != i) f_ij``
52 There is some freedom of choice in assigning atomic energies, i.e.
53 choosing a way to partition the total energy into atomic contributions.
55 We choose a symmetric approach (`bothways=True` in the neighbor list):
57 ``u_i = 1/2 sum_(j != i) u_ij``
59 The total energy of a system of atoms is then:
61 ``u = sum_i u_i = 1/2 sum_(i, j != i) u_ij``
63 Differentiating `u` with respect to `r_i` yields the force,
64 independent of the choice of partitioning.
66 ::
68 f_i = - d u / d r_i = - sum_ij d u_ij / d r_i
69 = - sum_ij d u_ij / d r_ij * d r_ij / d r_i
70 = sum_ij du_ij d_ij = sum_ij f_ij
72 This justifies calling `f_ij` pairwise forces.
74 The stress can be written as ( `(x)` denoting outer product):
76 ``sigma = 1/2 sum_(i, j != i) f_ij (x) d_ij = sum_i sigma_i ,``
77 with atomic contributions
79 ``sigma_i = 1/2 sum_(j != i) f_ij (x) d_ij``
81 Another consideration is the cutoff. We have to ensure that the potential
82 goes to zero smoothly as an atom moves across the cutoff threshold,
83 otherwise the potential is not continuous. In cases where the cutoff is
84 so large that u_ij is very small at the cutoff this is automatically
85 ensured, but in general, `u_ij(rc) != 0`.
87 This implementation offers two ways to deal with this:
89 Either, we shift the pairwise energy
91 ``u'_ij = u_ij - u_ij(rc)``
93 which ensures that it is precisely zero at the cutoff. However, this means
94 that the energy effectively depends on the cutoff, which might lead to
95 unexpected results! If this option is chosen, the forces discontinuously
96 jump to zero at the cutoff.
98 An alternative is to modify the pairwise potential by multiplying
99 it with a cutoff function that goes from 1 to 0 between an onset radius
100 ro and the cutoff rc. If the function is chosen suitably, it can also
101 smoothly push the forces down to zero, ensuring continuous forces as well.
102 In order for this to work well, the onset radius has to be set suitably,
103 typically around 2*sigma.
105 In this case, we introduce a modified pairwise potential:
107 ``u'_ij = fc * u_ij``
109 The pairwise forces have to be modified accordingly:
111 ``f'_ij = fc * f_ij + fc' * u_ij``
113 Where `fc' = d fc / d d_ij`.
115 This approach is taken from Jax-MD (https://github.com/google/jax-md),
116 which in turn is inspired by HOOMD Blue
117 (https://glotzerlab.engin.umich.edu/hoomd-blue/).
119 """
121 implemented_properties = ['energy', 'energies', 'forces', 'free_energy']
122 implemented_properties += ['stress', 'stresses'] # bulk properties
123 default_parameters = {
124 'epsilon': 1.0,
125 'sigma': 1.0,
126 'rc': None,
127 'ro': None,
128 'smooth': False,
129 }
130 nolabel = True
132 def __init__(self, **kwargs):
133 """
134 Parameters
135 ----------
136 sigma: float
137 The potential minimum is at 2**(1/6) * sigma, default 1.0
138 epsilon: float
139 The potential depth, default 1.0
140 rc: float, None
141 Cut-off for the NeighborList is set to 3 * sigma if None.
142 The energy is upshifted to be continuous at rc.
143 Default None
144 ro: float, None
145 Onset of cutoff function in 'smooth' mode. Defaults to 0.66 * rc.
146 smooth: bool, False
147 Cutoff mode. False means that the pairwise energy is simply shifted
148 to be 0 at r = rc, leading to the energy going to 0 continuously,
149 but the forces jumping to zero discontinuously at the cutoff.
150 True means that a smooth cutoff function is multiplied to the pairwise
151 energy that smoothly goes to 0 between ro and rc. Both energy and
152 forces are continuous in that case.
153 If smooth=True, make sure to check the tail of the
154 forces for kinks, ro might have to be adjusted to avoid distorting
155 the potential too much.
157 """
159 Calculator.__init__(self, **kwargs)
161 if self.parameters.rc is None:
162 self.parameters.rc = 3 * self.parameters.sigma
164 if self.parameters.ro is None:
165 self.parameters.ro = 0.66 * self.parameters.rc
167 self.nl = None
169 def calculate(
170 self,
171 atoms=None,
172 properties=None,
173 system_changes=all_changes,
174 ):
175 if properties is None:
176 properties = self.implemented_properties
178 Calculator.calculate(self, atoms, properties, system_changes)
180 natoms = len(self.atoms)
182 sigma = self.parameters.sigma
183 epsilon = self.parameters.epsilon
184 rc = self.parameters.rc
185 ro = self.parameters.ro
186 smooth = self.parameters.smooth
188 if self.nl is None or 'numbers' in system_changes:
189 self.nl = NeighborList(
190 [rc / 2] * natoms, self_interaction=False, bothways=True
191 )
193 self.nl.update(self.atoms)
195 positions = self.atoms.positions
196 cell = self.atoms.cell
198 # potential value at rc
199 e0 = 4 * epsilon * ((sigma / rc) ** 12 - (sigma / rc) ** 6)
201 energies = np.zeros(natoms)
202 forces = np.zeros((natoms, 3))
203 stresses = np.zeros((natoms, 3, 3))
205 for ii in range(natoms):
206 neighbors, offsets = self.nl.get_neighbors(ii)
207 cells = np.dot(offsets, cell)
209 # pointing *towards* neighbours
210 distance_vectors = positions[neighbors] + cells - positions[ii]
212 r2 = (distance_vectors ** 2).sum(1)
213 c6 = (sigma ** 2 / r2) ** 3
214 c6[r2 > rc ** 2] = 0.0
215 c12 = c6 ** 2
217 if smooth:
218 cutoff_fn = cutoff_function(r2, rc**2, ro**2)
219 d_cutoff_fn = d_cutoff_function(r2, rc**2, ro**2)
221 pairwise_energies = 4 * epsilon * (c12 - c6)
222 pairwise_forces = -24 * epsilon * (2 * c12 - c6) / r2 # du_ij
224 if smooth:
225 # order matters, otherwise the pairwise energy is already
226 # modified
227 pairwise_forces = (
228 cutoff_fn * pairwise_forces + 2 * d_cutoff_fn
229 * pairwise_energies
230 )
231 pairwise_energies *= cutoff_fn
232 else:
233 pairwise_energies -= e0 * (c6 != 0.0)
235 pairwise_forces = pairwise_forces[:, np.newaxis] * distance_vectors
237 energies[ii] += 0.5 * pairwise_energies.sum() # atomic energies
238 forces[ii] += pairwise_forces.sum(axis=0)
240 stresses[ii] += 0.5 * np.dot(
241 pairwise_forces.T, distance_vectors
242 ) # equivalent to outer product
244 # no lattice, no stress
245 if self.atoms.cell.rank == 3:
246 stresses = full_3x3_to_voigt_6_stress(stresses)
247 self.results['stress'] = stresses.sum(
248 axis=0) / self.atoms.get_volume()
249 self.results['stresses'] = stresses / self.atoms.get_volume()
251 energy = energies.sum()
252 self.results['energy'] = energy
253 self.results['energies'] = energies
255 self.results['free_energy'] = energy
257 self.results['forces'] = forces
260def cutoff_function(r, rc, ro):
261 """Smooth cutoff function.
263 Goes from 1 to 0 between ro and rc, ensuring
264 that u(r) = lj(r) * cutoff_function(r) is C^1.
266 Defined as 1 below ro, 0 above rc.
268 Note that r, rc, ro are all expected to be squared,
269 i.e. `r = r_ij^2`, etc.
271 Taken from https://github.com/google/jax-md.
273 """
275 return np.where(
276 r < ro,
277 1.0,
278 np.where(r < rc, (rc - r) ** 2 * (rc + 2 *
279 r - 3 * ro) / (rc - ro) ** 3, 0.0),
280 )
283def d_cutoff_function(r, rc, ro):
284 """Derivative of smooth cutoff function wrt r.
286 Note that `r = r_ij^2`, so for the derivative wrt to `r_ij`,
287 we need to multiply `2*r_ij`. This gives rise to the factor 2
288 above, the `r_ij` is cancelled out by the remaining derivative
289 `d r_ij / d d_ij`, i.e. going from scalar distance to distance vector.
290 """
292 return np.where(
293 r < ro,
294 0.0,
295 np.where(r < rc, 6 * (rc - r) * (ro - r) / (rc - ro) ** 3, 0.0),
296 )