 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) ( 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

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 2/3 * 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.

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 )