Coverage for hiphive/structure_generation/rattle.py: 85%
55 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-08-01 17:04 +0000
1"""
2Module for generating rattled structures. Rattle refers to displacing atoms
3with a normal distribution with zero mean and some standard deviation.
4"""
6import numpy as np
7from scipy.special import erf
8from ase import Atoms
9from ase.neighborlist import NeighborList
12def generate_rattled_structures(
13 atoms: Atoms,
14 n_structures: int,
15 rattle_std: float,
16 seed: int = 42,
17) -> list[Atoms]:
18 """Returns list of rattled configurations.
20 Displacements are drawn from normal distributions for each
21 Cartesian directions for each atom independently.
23 Warning
24 -------
25 Repeatedly calling this function *without* providing different
26 seeds will yield identical or correlated results. To avoid this
27 behavior it is recommended to specify a different seed for each
28 call to this function.
30 Parameters
31 ----------
32 atoms
33 Prototype structure.
34 n_structures
35 Number of structures to generate.
36 rattle_std
37 Rattle amplitude (standard deviation of the normal distribution).
38 seed
39 Seed for setting up NumPy random state from which random numbers are
40 generated.
42 Returns
43 -------
44 Generated structures.
45 """
46 rs = np.random.RandomState(seed)
47 N = len(atoms)
48 atoms_list = []
49 for _ in range(n_structures):
50 atoms_tmp = atoms.copy()
51 displacements = rs.normal(0.0, rattle_std, (N, 3))
52 atoms_tmp.positions += displacements
53 atoms_list.append(atoms_tmp)
54 return atoms_list
57def generate_mc_rattled_structures(
58 atoms: Atoms,
59 n_structures: int,
60 rattle_std: float,
61 d_min: float,
62 seed: int = 42,
63 **kwargs,
64) -> list[Atoms]:
65 r"""Returns list of Monte Carlo rattled configurations.
67 Rattling atom `i` is carried out as a Monte Carlo move that is
68 accepted with a probability determined from the minimum
69 interatomic distance :math:`d_{ij}`. If :math:`\min(d_{ij})` is
70 smaller than :math:`d_{min}` the move is only accepted with a low
71 probability.
73 This process is repeated for each atom a number of times meaning
74 the magnitude of the final displacements is not *directly*
75 connected to `rattle_std`.
77 Warning
78 -------
79 Repeatedly calling this function *without* providing different
80 seeds will yield identical or correlated results. To avoid this
81 behavior it is recommended to specify a different seed for each
82 call to this function.
84 Notes
85 ------
86 The procedure implemented here might not generate a symmetric
87 distribution for the displacements `kwargs` will be forwarded to
88 `mc_rattle` (see user guide for a detailed explanation)
90 The displacements generated will roughly be `n_iter**0.5 * rattle_std`
91 for small values of `n_iter`.
93 Parameters
94 ----------
95 atoms
96 Prototype structure.
97 n_structures
98 Number of structures to generate.
99 rattle_std
100 Rattle amplitude (standard deviation in normal distribution);
101 note this value is not connected to the final
102 average displacement for the structures.
103 d_min
104 Interatomic distance used for computing the probability for each rattle move.
105 seed
106 Seed for setting up NumPy random state from which random numbers are generated.
107 n_iter
108 Number of Monte Carlo cycles (iterations), larger number of iterations will
109 generate larger displacements (defaults to 10).
111 Returns
112 -------
113 Generated structures.
114 """
115 rs = np.random.RandomState(seed)
116 atoms_list = []
117 for _ in range(n_structures):
118 atoms_tmp = atoms.copy()
119 seed = rs.randint(1, 1000000000)
120 displacements = mc_rattle(atoms_tmp, rattle_std, d_min, seed=seed, **kwargs)
121 atoms_tmp.positions += displacements
122 atoms_list.append(atoms_tmp)
123 return atoms_list
126def _probability_mc_rattle(d: float, d_min: float, width: float):
127 """ Monte Carlo probability function as an error function.
129 Parameters
130 ----------
131 d
132 Value at which to evaluate function.
133 d_min
134 Center value for the error function.
135 width
136 Width of error function.
137 """
139 return (erf((d-d_min)/width) + 1.0) / 2
142def mc_rattle(
143 atoms: Atoms,
144 rattle_std: float,
145 d_min: float,
146 width: float = 0.1,
147 n_iter: int = 10,
148 max_attempts: int = 5000,
149 max_disp: float = 2.0,
150 active_atoms: list[int] = None,
151 nbr_cutoff: float = None,
152 seed: int = 42,
153) -> np.ndarray:
154 """Generate displacements using the Monte Carlo rattle method.
156 Parameters
157 ----------
158 atoms
159 Prototype structure.
160 rattle_std
161 Rattle amplitude (standard deviation in normal distribution).
162 d_min
163 Interatomic distance used for computing the probability for each rattle
164 move. Center position of the error function.
165 width
166 Width of the error function.
167 n_iter
168 Number of Monte Carlo cycle.
169 max_disp
170 Rattle moves that yields a displacement larger than max_disp will
171 always be rejected. This occurs rarely and is more used as a safety precaution
172 for not generating structures where two or more have swapped positions.
173 max_attempts
174 Limit for how many attempted rattle moves are allowed a single atom;
175 if this limit is reached an `Exception` is raised.
176 active_atoms
177 List of indices of atoms that should undergo Monte Carlo rattling.
178 nbr_cutoff
179 The cutoff used to construct the neighborlist used for checking
180 interatomic distances, defaults to `2 * d_min`.
181 seed
182 Seed for setting up NumPy random state from which random numbers are
183 generated.
185 Returns
186 -------
187 Atomic displacements (`Nx3`).
188 """
190 # setup
191 rs = np.random.RandomState(seed)
193 if nbr_cutoff is None: 193 ↛ 196line 193 didn't jump to line 196 because the condition on line 193 was always true
194 nbr_cutoff = 2 * d_min
196 if active_atoms is None: 196 ↛ 199line 196 didn't jump to line 199 because the condition on line 196 was always true
197 active_atoms = range(len(atoms))
199 atoms_rattle = atoms.copy()
200 reference_positions = atoms_rattle.get_positions()
201 nbr_list = NeighborList([nbr_cutoff/2]*len(atoms_rattle), skin=0.0,
202 self_interaction=False, bothways=True)
203 nbr_list.update(atoms_rattle)
205 # run Monte Carlo
206 for _ in range(n_iter):
207 for i in active_atoms:
208 i_nbrs = np.setdiff1d(nbr_list.get_neighbors(i)[0], [i])
209 for n in range(max_attempts): 209 ↛ 236line 209 didn't jump to line 236 because the loop on line 209 didn't complete
211 # generate displacement
212 delta_disp = rs.normal(0.0, rattle_std, 3)
213 atoms_rattle.positions[i] += delta_disp
214 disp_i = atoms_rattle.positions[i] - reference_positions[i]
216 # if total displacement of atom is greater than max_disp, then reject delta_disp
217 if np.linalg.norm(disp_i) > max_disp: 217 ↛ 219line 217 didn't jump to line 219 because the condition on line 217 was never true
218 # revert delta_disp
219 atoms_rattle[i].position -= delta_disp
220 continue
222 # compute min distance
223 if len(i_nbrs) == 0: 223 ↛ 224line 223 didn't jump to line 224 because the condition on line 223 was never true
224 min_distance = np.inf
225 else:
226 min_distance = np.min(atoms_rattle.get_distances(i, i_nbrs, mic=True))
228 # accept or reject delta_disp
229 if _probability_mc_rattle(min_distance, d_min, width) > rs.rand(): 229 ↛ 234line 229 didn't jump to line 234 because the condition on line 229 was always true
230 # accept delta_disp
231 break
232 else:
233 # revert delta_disp
234 atoms_rattle[i].position -= delta_disp
235 else:
236 raise Exception(f'Exceeded the maximum number of attempts for atom {i}')
237 displacements = atoms_rattle.positions - reference_positions
238 return displacements