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