Coverage for /builds/kinetik161/ase/ase/ga/soft_mutation.py: 90.20%
245 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Soft-mutation operator and associated tools"""
2import inspect
3import json
5import numpy as np
6from scipy.spatial.distance import cdist
8from ase.data import covalent_radii
9from ase.ga.offspring_creator import OffspringCreator
10from ase.ga.utilities import atoms_too_close, gather_atoms_by_tag
11from ase.neighborlist import NeighborList
14class TagFilter:
15 """Filter which constrains same-tag atoms to behave
16 like internally rigid moieties.
17 """
19 def __init__(self, atoms):
20 self.atoms = atoms
21 gather_atoms_by_tag(self.atoms)
22 self.tags = self.atoms.get_tags()
23 self.unique_tags = np.unique(self.tags)
24 self.n = len(self.unique_tags)
26 def get_positions(self):
27 all_pos = self.atoms.get_positions()
28 cop_pos = np.zeros((self.n, 3))
29 for i in range(self.n):
30 indices = np.where(self.tags == self.unique_tags[i])
31 cop_pos[i] = np.average(all_pos[indices], axis=0)
32 return cop_pos
34 def set_positions(self, positions, **kwargs):
35 cop_pos = self.get_positions()
36 all_pos = self.atoms.get_positions()
37 assert np.all(np.shape(positions) == np.shape(cop_pos))
38 for i in range(self.n):
39 indices = np.where(self.tags == self.unique_tags[i])
40 shift = positions[i] - cop_pos[i]
41 all_pos[indices] += shift
42 self.atoms.set_positions(all_pos, **kwargs)
44 def get_forces(self, *args, **kwargs):
45 f = self.atoms.get_forces()
46 forces = np.zeros((self.n, 3))
47 for i in range(self.n):
48 indices = np.where(self.tags == self.unique_tags[i])
49 forces[i] = np.sum(f[indices], axis=0)
50 return forces
52 def get_masses(self):
53 m = self.atoms.get_masses()
54 masses = np.zeros(self.n)
55 for i in range(self.n):
56 indices = np.where(self.tags == self.unique_tags[i])
57 masses[i] = np.sum(m[indices])
58 return masses
60 def __len__(self):
61 return self.n
64class PairwiseHarmonicPotential:
65 """Parent class for interatomic potentials of the type
66 E(r_ij) = 0.5 * k_ij * (r_ij - r0_ij) ** 2
67 """
69 def __init__(self, atoms, rcut=10.):
70 self.atoms = atoms
71 self.pos0 = atoms.get_positions()
72 self.rcut = rcut
74 # build neighborlist
75 nat = len(self.atoms)
76 self.nl = NeighborList([self.rcut / 2.] * nat, skin=0., bothways=True,
77 self_interaction=False)
78 self.nl.update(self.atoms)
80 self.calculate_force_constants()
82 def calculate_force_constants(self):
83 msg = 'Child class needs to define a calculate_force_constants() ' \
84 'method which computes the force constants and stores them ' \
85 'in self.force_constants (as a list which contains, for every ' \
86 'atom, a list of the atom\'s force constants with its neighbors.'
87 raise NotImplementedError(msg)
89 def get_forces(self, atoms):
90 pos = atoms.get_positions()
91 cell = atoms.get_cell()
92 forces = np.zeros_like(pos)
94 for i, p in enumerate(pos):
95 indices, offsets = self.nl.get_neighbors(i)
96 p = pos[indices] + np.dot(offsets, cell)
97 r = cdist(p, [pos[i]])
98 v = (p - pos[i]) / r
99 p0 = self.pos0[indices] + np.dot(offsets, cell)
100 r0 = cdist(p0, [self.pos0[i]])
101 dr = r - r0
102 forces[i] = np.dot(self.force_constants[i].T, dr * v)
104 return forces
107def get_number_of_valence_electrons(Z):
108 """Return the number of valence electrons for the element with
109 atomic number Z, simply based on its periodic table group.
110 """
111 groups = [[], [1, 3, 11, 19, 37, 55, 87], [2, 4, 12, 20, 38, 56, 88],
112 [21, 39, 57, 89]]
114 for i in range(9):
115 groups.append(i + np.array([22, 40, 72, 104]))
117 for i in range(6):
118 groups.append(i + np.array([5, 13, 31, 49, 81, 113]))
120 for i, group in enumerate(groups):
121 if Z in group:
122 nval = i if i < 13 else i - 10
123 break
124 else:
125 raise ValueError('Z=%d not included in this dataset.' % Z)
127 return nval
130class BondElectroNegativityModel(PairwiseHarmonicPotential):
131 """Pairwise harmonic potential where the force constants are
132 determined using the "bond electronegativity" model, see:
134 * `Lyakhov, Oganov, Valle, Comp. Phys. Comm. 181 (2010) 1623-1632`__
136 __ https://dx.doi.org/10.1016/j.cpc.2010.06.007
138 * `Lyakhov, Oganov, Phys. Rev. B 84 (2011) 092103`__
140 __ https://dx.doi.org/10.1103/PhysRevB.84.092103
141 """
143 def calculate_force_constants(self):
144 cell = self.atoms.get_cell()
145 pos = self.atoms.get_positions()
146 num = self.atoms.get_atomic_numbers()
147 nat = len(self.atoms)
148 nl = self.nl
150 # computing the force constants
151 s_norms = []
152 valence_states = []
153 r_cov = []
154 for i in range(nat):
155 indices, offsets = nl.get_neighbors(i)
156 p = pos[indices] + np.dot(offsets, cell)
157 r = cdist(p, [pos[i]])
158 r_ci = covalent_radii[num[i]]
159 s = 0.
160 for j, index in enumerate(indices):
161 d = r[j] - r_ci - covalent_radii[num[index]]
162 s += np.exp(-d / 0.37)
163 s_norms.append(s)
164 valence_states.append(get_number_of_valence_electrons(num[i]))
165 r_cov.append(r_ci)
167 self.force_constants = []
168 for i in range(nat):
169 indices, offsets = nl.get_neighbors(i)
170 p = pos[indices] + np.dot(offsets, cell)
171 r = cdist(p, [pos[i]])[:, 0]
172 fc = []
173 for j, ii in enumerate(indices):
174 d = r[j] - r_cov[i] - r_cov[ii]
175 chi_ik = 0.481 * valence_states[i] / (r_cov[i] + 0.5 * d)
176 chi_jk = 0.481 * valence_states[ii] / (r_cov[ii] + 0.5 * d)
177 cn_ik = s_norms[i] / np.exp(-d / 0.37)
178 cn_jk = s_norms[ii] / np.exp(-d / 0.37)
179 fc.append(np.sqrt(chi_ik * chi_jk / (cn_ik * cn_jk)))
180 self.force_constants.append(np.array(fc))
183class SoftMutation(OffspringCreator):
184 """Mutates the structure by displacing it along the lowest
185 (nonzero) frequency modes found by vibrational analysis, as in:
187 `Lyakhov, Oganov, Valle, Comp. Phys. Comm. 181 (2010) 1623-1632`__
189 __ https://dx.doi.org/10.1016/j.cpc.2010.06.007
191 As in the reference above, the next-lowest mode is used if the
192 structure has already been softmutated along the current-lowest
193 mode. This mutation hence acts in a deterministic way, in contrast
194 to most other genetic operators.
196 If you find this implementation useful in your work,
197 please consider citing:
199 `Van den Bossche, Gronbeck, Hammer, J. Chem. Theory Comput. 14 (2018)`__
201 __ https://dx.doi.org/10.1021/acs.jctc.8b00039
203 in addition to the paper mentioned above.
205 Parameters:
207 blmin: dict
208 The closest allowed interatomic distances on the form:
209 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
211 bounds: list
212 Lower and upper limits (in Angstrom) for the largest
213 atomic displacement in the structure. For a given mode,
214 the algorithm starts at zero amplitude and increases
215 it until either blmin is violated or the largest
216 displacement exceeds the provided upper bound).
217 If the largest displacement in the resulting structure
218 is lower than the provided lower bound, the mutant is
219 considered too similar to the parent and None is
220 returned.
222 calculator: ASE calculator object
223 The calculator to be used in the vibrational
224 analysis. The default is to use a calculator
225 based on pairwise harmonic potentials with force
226 constants from the "bond electronegativity"
227 model described in the reference above.
228 Any calculator with a working :func:`get_forces()`
229 method will work.
231 rcut: float
232 Cutoff radius in Angstrom for the pairwise harmonic
233 potentials.
235 used_modes_file: str or None
236 Name of json dump file where previously used
237 modes will be stored (and read). If None,
238 no such file will be used. Default is to use
239 the filename 'used_modes.json'.
241 use_tags: boolean
242 Whether to use the atomic tags to preserve molecular identity.
243 """
245 def __init__(self, blmin, bounds=[0.5, 2.0],
246 calculator=BondElectroNegativityModel, rcut=10.,
247 used_modes_file='used_modes.json', use_tags=False,
248 verbose=False):
249 OffspringCreator.__init__(self, verbose)
250 self.blmin = blmin
251 self.bounds = bounds
252 self.calc = calculator
253 self.rcut = rcut
254 self.used_modes_file = used_modes_file
255 self.use_tags = use_tags
256 self.descriptor = 'SoftMutation'
258 self.used_modes = {}
259 if self.used_modes_file is not None:
260 try:
261 self.read_used_modes(self.used_modes_file)
262 except OSError:
263 # file doesn't exist (yet)
264 pass
266 def _get_hessian(self, atoms, dx):
267 """Returns the Hessian matrix d2E/dxi/dxj using a first-order
268 central difference scheme with displacements dx.
269 """
270 N = len(atoms)
271 pos = atoms.get_positions()
272 hessian = np.zeros((3 * N, 3 * N))
274 for i in range(3 * N):
275 row = np.zeros(3 * N)
276 for direction in [-1, 1]:
277 disp = np.zeros(3)
278 disp[i % 3] = direction * dx
279 pos_disp = np.copy(pos)
280 pos_disp[i // 3] += disp
281 atoms.set_positions(pos_disp)
282 f = atoms.get_forces()
283 row += -1 * direction * f.flatten()
285 row /= (2. * dx)
286 hessian[i] = row
288 hessian += np.copy(hessian).T
289 hessian *= 0.5
290 atoms.set_positions(pos)
292 return hessian
294 def _calculate_normal_modes(self, atoms, dx=0.02, massweighing=False):
295 """Performs the vibrational analysis."""
296 hessian = self._get_hessian(atoms, dx)
297 if massweighing:
298 m = np.array([np.repeat(atoms.get_masses()**-0.5, 3)])
299 hessian *= (m * m.T)
301 eigvals, eigvecs = np.linalg.eigh(hessian)
302 modes = {eigval: eigvecs[:, i] for i, eigval in enumerate(eigvals)}
303 return modes
305 def animate_mode(self, atoms, mode, nim=30, amplitude=1.0):
306 """Returns an Atoms object showing an animation of the mode."""
307 pos = atoms.get_positions()
308 mode = mode.reshape(np.shape(pos))
309 animation = []
310 for i in range(nim):
311 newpos = pos + amplitude * mode * np.sin(i * 2 * np.pi / nim)
312 image = atoms.copy()
313 image.positions = newpos
314 animation.append(image)
315 return animation
317 def read_used_modes(self, filename):
318 """Read used modes from json file."""
319 with open(filename) as fd:
320 modes = json.load(fd)
321 self.used_modes = {int(k): modes[k] for k in modes}
322 return
324 def write_used_modes(self, filename):
325 """Dump used modes to json file."""
326 with open(filename, 'w') as fd:
327 json.dump(self.used_modes, fd)
328 return
330 def get_new_individual(self, parents):
331 f = parents[0]
333 indi = self.mutate(f)
334 if indi is None:
335 return indi, 'mutation: soft'
337 indi = self.initialize_individual(f, indi)
338 indi.info['data']['parents'] = [f.info['confid']]
340 return self.finalize_individual(indi), 'mutation: soft'
342 def mutate(self, atoms):
343 """Does the actual mutation."""
344 a = atoms.copy()
346 if inspect.isclass(self.calc):
347 assert issubclass(self.calc, PairwiseHarmonicPotential)
348 calc = self.calc(atoms, rcut=self.rcut)
349 else:
350 calc = self.calc
351 a.calc = calc
353 if self.use_tags:
354 a = TagFilter(a)
356 pos = a.get_positions()
357 modes = self._calculate_normal_modes(a)
359 # Select the mode along which we want to move the atoms;
360 # The first 3 translational modes as well as previously
361 # applied modes are discarded.
363 keys = np.array(sorted(modes))
364 index = 3
365 confid = atoms.info['confid']
366 if confid in self.used_modes:
367 while index in self.used_modes[confid]:
368 index += 1
369 self.used_modes[confid].append(index)
370 else:
371 self.used_modes[confid] = [index]
373 if self.used_modes_file is not None:
374 self.write_used_modes(self.used_modes_file)
376 key = keys[index]
377 mode = modes[key].reshape(np.shape(pos))
379 # Find a suitable amplitude for translation along the mode;
380 # at every trial amplitude both positive and negative
381 # directions are tried.
383 mutant = atoms.copy()
384 amplitude = 0.
385 increment = 0.1
386 direction = 1
387 largest_norm = np.max(np.apply_along_axis(np.linalg.norm, 1, mode))
389 def expand(atoms, positions):
390 if isinstance(atoms, TagFilter):
391 a.set_positions(positions)
392 return a.atoms.get_positions()
393 else:
394 return positions
396 while amplitude * largest_norm < self.bounds[1]:
397 pos_new = pos + direction * amplitude * mode
398 pos_new = expand(a, pos_new)
399 mutant.set_positions(pos_new)
400 mutant.wrap()
401 too_close = atoms_too_close(mutant, self.blmin,
402 use_tags=self.use_tags)
403 if too_close:
404 amplitude -= increment
405 pos_new = pos + direction * amplitude * mode
406 pos_new = expand(a, pos_new)
407 mutant.set_positions(pos_new)
408 mutant.wrap()
409 break
411 if direction == 1:
412 direction = -1
413 else:
414 direction = 1
415 amplitude += increment
417 if amplitude * largest_norm < self.bounds[0]:
418 mutant = None
420 return mutant