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

1# Contains the get_clusters function which generates clusters 

2 

3import numpy as np 

4import itertools 

5from collections import defaultdict 

6 

7from .utilities import BiMap 

8from ..input_output.logging_tools import logger 

9 

10logger = logger.getChild('get_clusters') 

11 

12 

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

14def get_clusters(atoms, cutoffs, nPrim, multiplicity=True, 

15 use_geometrical_order=False): 

16 """Generate a list of all clusters in the atoms object which includes the 

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

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

19 

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

21 be generated. This is useful when doing force constants but not so much for 

22 cluster expansions. 

23 

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

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

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

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

28 many body interactions decrease fast with cutoff but anharmonic 

29 interactions can be quite long ranged. 

30 

31 Parameters 

32 ---------- 

33 atoms : ase.Atoms 

34 can be a general atoms object but must have pbc=False. 

35 cutoffs : dict 

36 the keys specify the order while the values specify the cutoff radii 

37 multiplicity : bool 

38 includes clusters where same atom appears more than once 

39 geometrical_order : bool 

40 specifies if the geometrical order should be used as cutoff_order, 

41 otherwise the normal order of the cluster is used 

42 

43 Returns 

44 ------- 

45 list(tuple(int)) 

46 a list of clusters where each entry is a tuple of indices, 

47 which refer to the atoms in the input supercell 

48 """ 

49 

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

51 cluster_dict = defaultdict(list) 

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

53 for i in range(nPrim): 

54 for order in cutoffs.orders: 

55 cluster = (i,) * order 

56 cluster_dict[order].append(cluster) 

57 

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

59 for nbody in cutoffs.nbodies: 

60 cutoff = cutoffs.max_nbody_cutoff(nbody) 

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

62 nbody_clusters, nbody_cutoffs = generate_geometrical_clusters(atoms, nPrim, cutoff, nbody) 

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

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

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

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

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

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

69 

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

71 cluster_list = BiMap() 

72 for key in sorted(cluster_dict): 

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

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

75 cluster_list.append(cluster) 

76 return cluster_list 

77 

78 

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

80 neighbor_matrix, distance_matrix = create_neighbor_matrices(atoms, cutoff) 

81 clusters, cutoffs = [], [] 

82 i, j = 0, 0 

83 # The clusters are generated in lexicographical order 

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

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

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

87 # we can break 

88 if cluster[0] >= n_prim: 

89 break 

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

91 # things up 

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

93 continue 

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

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

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

97 break 

98 else: 

99 clusters.append(cluster) 

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

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

102 return clusters, cutoffs 

103 

104 

105def create_neighbor_matrices(atoms, cutoff): 

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

107 neighbor_matrix = distance_matrix < cutoff 

108 return neighbor_matrix, distance_matrix 

109 

110 

111def extend_cluster(cluster, order): 

112 clusters = [] 

113 cluster = tuple(cluster) 

114 nbody = len(cluster) 

115 r = order - nbody 

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

117 new_cluster = sorted(cluster + tup) 

118 clusters.append(tuple(new_cluster)) 

119 return clusters