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