Coverage for hiphive/core/clusters.py: 100%

1# Contains the get_clusters function which generates clusters 

2 

3import itertools 

4from collections import defaultdict 

5 

6import ase 

7import numpy as np 

8from .utilities import BiMap 

9from ..input_output.logging_tools import logger 

10 

11logger = logger.getChild('get_clusters') 

12 

13 

14# TODO: This function could be made a bit more general 

15def get_clusters( 

16 atoms: ase.Atoms, 

17 cutoffs: dict[float], 

18 n_prim: int, 

19 multiplicity: bool = True, 

20 use_geometrical_order: bool = False, 

21) -> list[tuple[int]]: 

22 """Generates a list of all clusters in the atoms object which includes the 

23 center atoms with positions within the cell metric. The cutoff determines 

24 up to which order and range clusters should be generated. 

25 

26 With multiplicity set to True clusters like `[0, 0]` and `[3, 3, 4]` etc will 

27 be generated. This is useful when working with force constants but not so much 

28 for cluster expansions. 

29 

30 The geometrical order is the total number of different atoms in the 

31 cluster. `[0, 0, 1]` would have geometrical order 2 and `[1, 2, 3, 4]` would 

32 have order 4. If the key word is True the cutoff criteria will be based on 

33 the geometrical order of the cluster. This is based on the observation that 

34 many body interactions decrease fast with cutoff but anharmonic 

35 interactions can be quite long ranged. 

36 

37 Parameters 

38 ---------- 

39 atoms 

40 Atomic structure; can be a general ASE:class:`Atoms` object but must have `pbc=False`. 

41 cutoffs 

42 Cutoffs by order; the keys specify the order while the values specify the cutoff radii. 

43 n_prim 

44 Number of atoms in the primitive cell. 

45 multiplicity 

46 Include clusters in which same atom appears more than once. 

47 use_geometrical_order 

48 Specify if the geometrical order should be used as cutoff order, 

49 otherwise the normal order of the cluster is used. 

50 

51 Returns 

52 ------- 

53 A list of clusters where each entry is a tuple of indices, 

54 which refer to the atoms in the input supercell. 

55 """ 

56 

57 logger.debug('Generating clusters...') 

58 cluster_dict = defaultdict(list) 

59 # Generate all on-site clusters of all orders (1-body) 

60 for i in range(n_prim): 

61 for order in cutoffs.orders: 

62 cluster = (i,) * order 

63 cluster_dict[order].append(cluster) 

64 

65 # Generate all 2-body clusters and above in order 

66 for nbody in cutoffs.nbodies: 

67 cutoff = cutoffs.max_nbody_cutoff(nbody) 

68 # Generate all n-body, order n clusters compatible with the cutoff 

69 nbody_clusters, nbody_cutoffs = generate_geometrical_clusters(atoms, n_prim, cutoff, nbody) 

70 for order in range(nbody, cutoffs.max_nbody_order(nbody) + 1): 

71 for cluster, cutoff in zip(nbody_clusters, nbody_cutoffs): 

72 # If the cutoff of the n-body cluster is compatible with order (order > n) then 

73 # extend the n-body cluster to higher order (e.g. nbody=3, order=6: ijk -> iijkkk) 

74 if cutoff < cutoffs.get_cutoff(nbody=nbody, order=order): 

75 cluster_dict[order].extend(extend_cluster(cluster, order)) 

76 

77 # The clusters are saved in a BiMap structure which allows for fast lookups 

78 cluster_list = BiMap() 

79 for key in sorted(cluster_dict): 

80 # For each order the clusters are saved in lexicographical order 

81 for cluster in sorted(cluster_dict[key]): 

82 cluster_list.append(cluster) 

83 return cluster_list 

84 

85 

86def generate_geometrical_clusters(atoms, n_prim, cutoff, order): 

87 neighbor_matrix, distance_matrix = create_neighbor_matrices(atoms, cutoff) 

88 clusters, cutoffs = [], [] 

89 i, j = 0, 0 

90 # The clusters are generated in lexicographical order 

91 for cluster in itertools.combinations(range(len(atoms)), r=order): 

92 # If the first atom in the cluster has an index higher or equal to the number of atoms in 

93 # the primitive cell then no upcoming cluster will have an atom in the primitive cell, thus 

94 # we can break 

95 if cluster[0] >= n_prim: 

96 break 

97 # if the last cluster failed on index i, j we start by checking this index again to speed 

98 # things up 

99 if not neighbor_matrix[cluster[i], cluster[j]]: 

100 continue 

101 # loop through all pairs in the cluster and check so that they are neighbors 

102 for i, j in itertools.combinations(range(order), r=2): 

103 if not neighbor_matrix[cluster[i], cluster[j]]: 

104 break 

105 else: 

106 clusters.append(cluster) 

107 # We also note the cutoff each cluster is compatible with 

108 cutoffs.append(np.max(distance_matrix[cluster, :][:, cluster])) 

109 return clusters, cutoffs 

110 

111 

112def create_neighbor_matrices(atoms, cutoff): 

113 distance_matrix = atoms.get_all_distances(mic=False) # or True? 

114 neighbor_matrix = distance_matrix < cutoff 

115 return neighbor_matrix, distance_matrix 

116 

117 

118def extend_cluster(cluster, order): 

119 clusters = [] 

120 cluster = tuple(cluster) 

121 nbody = len(cluster) 

122 r = order - nbody 

123 for tup in itertools.combinations_with_replacement(cluster, r): 

124 new_cluster = sorted(cluster + tup) 

125 clusters.append(tuple(new_cluster)) 

126 return clusters