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