Calculates a python reference value for the Coulomb matrix.
(system, n_atoms_max, permutation, sigma=None)
| 25 | # ============================================================================= |
| 26 | # Utilities |
| 27 | def cm_python(system, n_atoms_max, permutation, sigma=None): |
| 28 | """Calculates a python reference value for the Coulomb matrix.""" |
| 29 | pos = system.get_positions() |
| 30 | n = len(system) |
| 31 | distances = np.linalg.norm(pos[:, None, :] - pos[None, :, :], axis=-1) |
| 32 | q = system.get_atomic_numbers() |
| 33 | qiqj = q[None, :] * q[:, None] |
| 34 | np.fill_diagonal(distances, 1) |
| 35 | cm = qiqj / distances |
| 36 | np.fill_diagonal(cm, 0.5 * q**2.4) |
| 37 | random_state = RandomState(42) |
| 38 | |
| 39 | # Permutation option |
| 40 | if permutation == "eigenspectrum": |
| 41 | eigenvalues = np.linalg.eigvalsh(cm) |
| 42 | abs_values = np.absolute(eigenvalues) |
| 43 | sorted_indices = np.argsort(abs_values)[::-1] |
| 44 | eigenvalues = eigenvalues[sorted_indices] |
| 45 | padded = np.zeros((n_atoms_max)) |
| 46 | padded[:n] = eigenvalues |
| 47 | else: |
| 48 | if permutation == "sorted_l2": |
| 49 | norms = np.linalg.norm(cm, axis=1) |
| 50 | sorted_indices = np.argsort(norms, axis=0)[::-1] |
| 51 | cm = cm[sorted_indices] |
| 52 | cm = cm[:, sorted_indices] |
| 53 | elif permutation == "random": |
| 54 | norms = np.linalg.norm(cm, axis=1) |
| 55 | noise_norm_vector = random_state.normal(norms, sigma) |
| 56 | indexlist = np.argsort(noise_norm_vector) |
| 57 | indexlist = indexlist[::-1] # Order highest to lowest |
| 58 | cm = cm[indexlist][:, indexlist] |
| 59 | elif permutation == "none": |
| 60 | pass |
| 61 | else: |
| 62 | raise ValueError("Unknown permutation option") |
| 63 | |
| 64 | # Flattening |
| 65 | padded = np.zeros((n_atoms_max, n_atoms_max)) |
| 66 | padded[:n, :n] = cm |
| 67 | padded = padded.flatten() |
| 68 | |
| 69 | return padded |
| 70 | |
| 71 | |
| 72 | def coulomb_matrix(**kwargs): |
no test coverage detected