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