Coverage for /builds/debichem-team/python-ase/ase/calculators/lj.py : 100.00%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import numpy as np
3from ase.neighborlist import NeighborList
4from ase.calculators.calculator import Calculator, all_changes
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) ( 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, indepedent of the
62 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), which in
114 turn is inspired by HOOMD Blue (https://glotzerlab.engin.umich.edu/hoomd-blue/).
116 """
118 implemented_properties = ['energy', 'energies', 'forces', 'free_energy']
119 implemented_properties += ['stress', 'stresses'] # bulk properties
120 default_parameters = {
121 'epsilon': 1.0,
122 'sigma': 1.0,
123 'rc': None,
124 'ro': None,
125 'smooth': False,
126 }
127 nolabel = True
129 def __init__(self, **kwargs):
130 """
131 Parameters
132 ----------
133 sigma: float
134 The potential minimum is at 2**(1/6) * sigma, default 1.0
135 epsilon: float
136 The potential depth, default 1.0
137 rc: float, None
138 Cut-off for the NeighborList is set to 3 * sigma if None.
139 The energy is upshifted to be continuous at rc.
140 Default None
141 ro: float, None
142 Onset of cutoff function in 'smooth' mode. Defaults to 2/3 * rc.
143 smooth: bool, False
144 Cutoff mode. False means that the pairwise energy is simply shifted
145 to be 0 at r = rc, leading to the energy going to 0 continuously,
146 but the forces jumping to zero discontinuously at the cutoff.
147 True means that a smooth cutoff function is multiplied to the pairwise
148 energy that smoothly goes to 0 between ro and rc. Both energy and
149 forces are continuous in that case.
150 If smooth=True, make sure to check the tail of the forces for kinks, ro
151 might have to be adjusted to avoid distorting the potential too much.
153 """
155 Calculator.__init__(self, **kwargs)
157 if self.parameters.rc is None:
158 self.parameters.rc = 3 * self.parameters.sigma
160 if self.parameters.ro is None:
161 self.parameters.ro = 0.66 * self.parameters.rc
163 self.nl = None
165 def calculate(
166 self,
167 atoms=None,
168 properties=None,
169 system_changes=all_changes,
170 ):
171 if properties is None:
172 properties = self.implemented_properties
174 Calculator.calculate(self, atoms, properties, system_changes)
176 natoms = len(self.atoms)
178 sigma = self.parameters.sigma
179 epsilon = self.parameters.epsilon
180 rc = self.parameters.rc
181 ro = self.parameters.ro
182 smooth = self.parameters.smooth
184 if self.nl is None or 'numbers' in system_changes:
185 self.nl = NeighborList(
186 [rc / 2] * natoms, self_interaction=False, bothways=True
187 )
189 self.nl.update(self.atoms)
191 positions = self.atoms.positions
192 cell = self.atoms.cell
194 # potential value at rc
195 e0 = 4 * epsilon * ((sigma / rc) ** 12 - (sigma / rc) ** 6)
197 energies = np.zeros(natoms)
198 forces = np.zeros((natoms, 3))
199 stresses = np.zeros((natoms, 3, 3))
201 for ii in range(natoms):
202 neighbors, offsets = self.nl.get_neighbors(ii)
203 cells = np.dot(offsets, cell)
205 # pointing *towards* neighbours
206 distance_vectors = positions[neighbors] + cells - positions[ii]
208 r2 = (distance_vectors ** 2).sum(1)
209 c6 = (sigma ** 2 / r2) ** 3
210 c6[r2 > rc ** 2] = 0.0
211 c12 = c6 ** 2
213 if smooth:
214 cutoff_fn = cutoff_function(r2, rc**2, ro**2)
215 d_cutoff_fn = d_cutoff_function(r2, rc**2, ro**2)
217 pairwise_energies = 4 * epsilon * (c12 - c6)
218 pairwise_forces = -24 * epsilon * (2 * c12 - c6) / r2 # du_ij
220 if smooth:
221 # order matters, otherwise the pairwise energy is already modified
222 pairwise_forces = (
223 cutoff_fn * pairwise_forces + 2 * d_cutoff_fn * pairwise_energies
224 )
225 pairwise_energies *= cutoff_fn
226 else:
227 pairwise_energies -= e0 * (c6 != 0.0)
229 pairwise_forces = pairwise_forces[:, np.newaxis] * distance_vectors
231 energies[ii] += 0.5 * pairwise_energies.sum() # atomic energies
232 forces[ii] += pairwise_forces.sum(axis=0)
234 stresses[ii] += 0.5 * np.dot(
235 pairwise_forces.T, distance_vectors
236 ) # equivalent to outer product
238 # no lattice, no stress
239 if self.atoms.cell.rank == 3:
240 stresses = full_3x3_to_voigt_6_stress(stresses)
241 self.results['stress'] = stresses.sum(axis=0) / self.atoms.get_volume()
242 self.results['stresses'] = stresses / self.atoms.get_volume()
244 energy = energies.sum()
245 self.results['energy'] = energy
246 self.results['energies'] = energies
248 self.results['free_energy'] = energy
250 self.results['forces'] = forces
253def cutoff_function(r, rc, ro):
254 """Smooth cutoff function.
256 Goes from 1 to 0 between ro and rc, ensuring
257 that u(r) = lj(r) * cutoff_function(r) is C^1.
259 Defined as 1 below ro, 0 above rc.
261 Note that r, rc, ro are all expected to be squared,
262 i.e. `r = r_ij^2`, etc.
264 Taken from https://github.com/google/jax-md.
266 """
268 return np.where(
269 r < ro,
270 1.0,
271 np.where(r < rc, (rc - r) ** 2 * (rc + 2 * r - 3 * ro) / (rc - ro) ** 3, 0.0),
272 )
275def d_cutoff_function(r, rc, ro):
276 """Derivative of smooth cutoff function wrt r.
278 Note that `r = r_ij^2`, so for the derivative wrt to `r_ij`,
279 we need to multiply `2*r_ij`. This gives rise to the factor 2
280 above, the `r_ij` is cancelled out by the remaining derivative
281 `d r_ij / d d_ij`, i.e. going from scalar distance to distance vector.
282 """
284 return np.where(
285 r < ro,
286 0.0,
287 np.where(r < rc, 6 * (rc - r) * (ro - r) / (rc - ro) ** 3, 0.0),
288 )