Coverage for /builds/kinetik161/ase/ase/build/supercells.py: 49.51%
103 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Helper functions for creating supercells."""
3import numpy as np
4from ase import Atoms
7class SupercellError(Exception):
8 """Use if construction of supercell fails"""
11def get_deviation_from_optimal_cell_shape(cell, target_shape="sc", norm=None):
12 r"""
13 Calculates the deviation of the given cell metric from the ideal
14 cell metric defining a certain shape. Specifically, the function
15 evaluates the expression `\Delta = || Q \mathbf{h} -
16 \mathbf{h}_{target}||_2`, where `\mathbf{h}` is the input
17 metric (*cell*) and `Q` is a normalization factor (*norm*)
18 while the target metric `\mathbf{h}_{target}` (via
19 *target_shape*) represent simple cubic ('sc') or face-centered
20 cubic ('fcc') cell shapes.
22 Parameters:
24 cell: 2D array of floats
25 Metric given as a (3x3 matrix) of the input structure.
26 target_shape: str
27 Desired supercell shape. Can be 'sc' for simple cubic or
28 'fcc' for face-centered cubic.
29 norm: float
30 Specify the normalization factor. This is useful to avoid
31 recomputing the normalization factor when computing the
32 deviation for a series of P matrices.
34 """
36 if target_shape in ["sc", "simple-cubic"]:
37 target_metric = np.eye(3)
38 elif target_shape in ["fcc", "face-centered cubic"]:
39 target_metric = 0.5 * np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
40 if not norm:
41 norm = (np.linalg.det(cell) / np.linalg.det(target_metric))**(-1.0 / 3)
42 return np.linalg.norm(norm * cell - target_metric)
45def find_optimal_cell_shape(
46 cell,
47 target_size,
48 target_shape,
49 lower_limit=-2,
50 upper_limit=2,
51 verbose=False,
52):
53 """Returns the transformation matrix that produces a supercell
54 corresponding to *target_size* unit cells with metric *cell* that
55 most closely approximates the shape defined by *target_shape*.
57 Parameters:
59 cell: 2D array of floats
60 Metric given as a (3x3 matrix) of the input structure.
61 target_size: integer
62 Size of desired super cell in number of unit cells.
63 target_shape: str
64 Desired supercell shape. Can be 'sc' for simple cubic or
65 'fcc' for face-centered cubic.
66 lower_limit: int
67 Lower limit of search range.
68 upper_limit: int
69 Upper limit of search range.
70 verbose: bool
71 Set to True to obtain additional information regarding
72 construction of transformation matrix.
74 """
76 # Set up target metric
77 if target_shape in ["sc", "simple-cubic"]:
78 target_metric = np.eye(3)
79 elif target_shape in ["fcc", "face-centered cubic"]:
80 target_metric = 0.5 * np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]],
81 dtype=float)
82 if verbose:
83 print("target metric (h_target):")
84 print(target_metric)
86 # Normalize cell metric to reduce computation time during looping
87 norm = (target_size * np.linalg.det(cell) /
88 np.linalg.det(target_metric))**(-1.0 / 3)
89 norm_cell = norm * cell
90 if verbose:
91 print("normalization factor (Q): %g" % norm)
93 # Approximate initial P matrix
94 ideal_P = np.dot(target_metric, np.linalg.inv(norm_cell))
95 if verbose:
96 print("idealized transformation matrix:")
97 print(ideal_P)
98 starting_P = np.array(np.around(ideal_P, 0), dtype=int)
99 if verbose:
100 print("closest integer transformation matrix (P_0):")
101 print(starting_P)
103 # Prepare run.
104 from itertools import product
106 best_score = 1e6
107 optimal_P = None
108 for dP in product(range(lower_limit, upper_limit + 1), repeat=9):
109 dP = np.array(dP, dtype=int).reshape(3, 3)
110 P = starting_P + dP
111 if int(np.around(np.linalg.det(P), 0)) != target_size:
112 continue
113 score = get_deviation_from_optimal_cell_shape(
114 np.dot(P, norm_cell), target_shape=target_shape, norm=1.0)
115 if score < best_score:
116 best_score = score
117 optimal_P = P
119 if optimal_P is None:
120 print("Failed to find a transformation matrix.")
121 return None
123 # Finalize.
124 if verbose:
125 print("smallest score (|Q P h_p - h_target|_2): %f" % best_score)
126 print("optimal transformation matrix (P_opt):")
127 print(optimal_P)
128 print("supercell metric:")
129 print(np.round(np.dot(optimal_P, cell), 4))
130 print("determinant of optimal transformation matrix: %g" %
131 np.linalg.det(optimal_P))
132 return optimal_P
135def make_supercell(prim, P, *, wrap=True, order="cell-major", tol=1e-5):
136 r"""Generate a supercell by applying a general transformation (*P*) to
137 the input configuration (*prim*).
139 The transformation is described by a 3x3 integer matrix
140 `\mathbf{P}`. Specifically, the new cell metric
141 `\mathbf{h}` is given in terms of the metric of the input
142 configuration `\mathbf{h}_p` by `\mathbf{P h}_p =
143 \mathbf{h}`.
145 Parameters:
147 prim: ASE Atoms object
148 Input configuration.
149 P: 3x3 integer matrix
150 Transformation matrix `\mathbf{P}`.
151 wrap: bool
152 wrap in the end
153 order: str (default: "cell-major")
154 how to order the atoms in the supercell
156 "cell-major":
157 [atom1_shift1, atom2_shift1, ..., atom1_shift2, atom2_shift2, ...]
158 i.e. run first over all the atoms in cell1 and then move to cell2.
160 "atom-major":
161 [atom1_shift1, atom1_shift2, ..., atom2_shift1, atom2_shift2, ...]
162 i.e. run first over atom1 in all the cells and then move to atom2.
163 This may be the order preferred by most VASP users.
165 tol: float
166 tolerance for wrapping
167 """
169 supercell_matrix = P
170 supercell = clean_matrix(supercell_matrix @ prim.cell)
172 # cartesian lattice points
173 lattice_points_frac = lattice_points_in_supercell(supercell_matrix)
174 lattice_points = np.dot(lattice_points_frac, supercell)
175 N = len(lattice_points)
177 if order == "cell-major":
178 shifted = prim.positions[None, :, :] + lattice_points[:, None, :]
179 elif order == "atom-major":
180 shifted = prim.positions[:, None, :] + lattice_points[None, :, :]
181 else:
182 raise ValueError(f"invalid order: {order}")
183 shifted_reshaped = shifted.reshape(-1, 3)
185 superatoms = Atoms(positions=shifted_reshaped,
186 cell=supercell,
187 pbc=prim.pbc)
189 # Copy over any other possible arrays, inspired by atoms.__imul__
190 for name, arr in prim.arrays.items():
191 if name == "positions":
192 # This was added during construction of the super cell
193 continue
194 shape = (N * arr.shape[0], *arr.shape[1:])
195 if order == "cell-major":
196 new_arr = np.repeat(arr[None, :], N, axis=0).reshape(shape)
197 elif order == "atom-major":
198 new_arr = np.repeat(arr[:, None], N, axis=1).reshape(shape)
199 superatoms.set_array(name, new_arr)
201 # check number of atoms is correct
202 n_target = abs(int(np.round(np.linalg.det(supercell_matrix) * len(prim))))
203 if n_target != len(superatoms):
204 msg = "Number of atoms in supercell: {}, expected: {}".format(
205 n_target, len(superatoms))
206 raise SupercellError(msg)
208 if wrap:
209 superatoms.wrap(eps=tol)
211 return superatoms
214def lattice_points_in_supercell(supercell_matrix):
215 """Find all lattice points contained in a supercell.
217 Adapted from pymatgen, which is available under MIT license:
218 The MIT License (MIT) Copyright (c) 2011-2012 MIT & The Regents of the
219 University of California, through Lawrence Berkeley National Laboratory
220 """
222 diagonals = np.array([
223 [0, 0, 0],
224 [0, 0, 1],
225 [0, 1, 0],
226 [0, 1, 1],
227 [1, 0, 0],
228 [1, 0, 1],
229 [1, 1, 0],
230 [1, 1, 1],
231 ])
232 d_points = np.dot(diagonals, supercell_matrix)
234 mins = np.min(d_points, axis=0)
235 maxes = np.max(d_points, axis=0) + 1
237 ar = np.arange(mins[0], maxes[0])[:, None] * np.array([1, 0, 0])[None, :]
238 br = np.arange(mins[1], maxes[1])[:, None] * np.array([0, 1, 0])[None, :]
239 cr = np.arange(mins[2], maxes[2])[:, None] * np.array([0, 0, 1])[None, :]
241 all_points = ar[:, None, None] + br[None, :, None] + cr[None, None, :]
242 all_points = all_points.reshape((-1, 3))
244 frac_points = np.dot(all_points, np.linalg.inv(supercell_matrix))
246 tvects = frac_points[np.all(frac_points < 1 - 1e-10, axis=1)
247 & np.all(frac_points >= -1e-10, axis=1)]
248 assert len(tvects) == round(abs(np.linalg.det(supercell_matrix)))
249 return tvects
252def clean_matrix(matrix, eps=1e-12):
253 """ clean from small values"""
254 matrix = np.array(matrix)
255 for ij in np.ndindex(matrix.shape):
256 if abs(matrix[ij]) < eps:
257 matrix[ij] = 0
258 return matrix