Coverage for mchammer/observers/site_occupancy_observer.py: 100%

1from ase import Atoms 

2from icet import ClusterSpace 

3from icet.core.structure import Structure 

4from mchammer.observers.base_observer import BaseObserver 

5from typing import List, Dict 

6import numpy as np 

7 

8 

9class SiteOccupancyObserver(BaseObserver): 

10 """ 

11 This class represents a site occupation factor (SOF) observer. 

12 

13 A SOF observer allows to compute the site occupation factors along the 

14 trajectory sampled by a Monte Carlo (MC) simulation. 

15 

16 Parameters 

17 ---------- 

18 cluster_space 

19 Cluster space from which the allowed species are extracted. 

20 structure 

21 Supercell consistent with primitive structure in :attr:`cluster_space`. 

22 Used to determine which species are allowed on each site. 

23 sites 

24 Dictionary containing lists of sites that are to be considered. 

25 The keys will be taken as the names of the sites; the indices refer 

26 to the primitive structure associated with the cluster space. 

27 interval 

28 Observation interval. Defaults to ``None`` meaning that if the 

29 observer is used in a Monte Carlo simulation, then the 

30 :class:`Ensemble` object will determine the interval. 

31 

32 Example 

33 ------- 

34 The following snippet illustrates how to use the site occupancy factor (SOF) 

35 observer in a Monte Carlo simulation of a surface slab. Here, the SOF 

36 observer is used to monitor the concentrations of different species at the 

37 surface, the first subsurface layer, and the remaining 'bulk'. A minimal 

38 cluster expansion is used with slightly modified surface interactions in 

39 order to obtain an example that can be run without much ado. In practice, 

40 one should of course use a proper cluster expansion:: 

41 

42 >>> from ase.build import fcc111 

43 >>> from icet import ClusterExpansion, ClusterSpace 

44 >>> from mchammer.calculators import ClusterExpansionCalculator 

45 >>> from mchammer.ensembles import CanonicalEnsemble 

46 >>> from mchammer.observers import SiteOccupancyObserver 

47 

48 >>> # prepare reference structure 

49 >>> prim = fcc111('Au', size=(1, 1, 10), vacuum=10.0) 

50 >>> prim.translate((0.1, 0.1, 0.0)) 

51 >>> prim.wrap() 

52 >>> prim.pbc = True # icet requires pbc in all directions 

53 

54 >>> # prepare cluster expansion 

55 >>> cs = ClusterSpace(prim, cutoffs=[3.7], chemical_symbols=['Ag', 'Au']) 

56 >>> params = [0] + 5 * [0] + 10 * [0.1] 

57 >>> params[1] = 0.01 

58 >>> params[6] = 0.12 

59 >>> ce = ClusterExpansion(cs, params) 

60 >>> print(ce) 

61 

62 >>> # prepare initial configuration based on a 2x2 supercell 

63 >>> structure = prim.repeat((2, 2, 1)) 

64 >>> for k in range(20): 

65 >>> structure[k].symbol = 'Ag' 

66 

67 >>> # set up MC simulation 

68 >>> calc = ClusterExpansionCalculator(structure, ce) 

69 >>> mc = CanonicalEnsemble(structure=structure, calculator=calc, temperature=600, 

70 ... dc_filename='myrun_sof.dc') 

71 

72 >>> # set up observer and attach it to the MC simulation 

73 >>> sites = {'surface': [0, 9], 'subsurface': [1, 8], 

74 ... 'bulk': list(range(2, 8))} 

75 >>> sof = SiteOccupancyObserver(cs, structure, sites, interval=len(structure)) 

76 >>> mc.attach_observer(sof) 

77 

78 >>> # run 1000 trial steps 

79 >>> mc.run(1000) 

80 

81 After having run this snippet one can access the SOFs via the data 

82 container:: 

83 

84 >>> print(mc.data_container.data) 

85 """ 

86 

87 def __init__(self, cluster_space: ClusterSpace, 

88 structure: Atoms, 

89 sites: Dict[str, List[int]], 

90 interval: int = None) -> None: 

91 super().__init__(interval=interval, return_type=dict, tag='SiteOccupancyObserver') 

92 

93 self._sites = {site: sorted(indices) 

94 for site, indices in sites.items()} 

95 

96 self._set_allowed_species(cluster_space, structure) 

97 

98 def _set_allowed_species(self, 

99 cluster_space: ClusterSpace, 

100 structure: Atoms): 

101 """ 

102 Set the allowed species for the selected sites in the Atoms object 

103 

104 Parameters 

105 ---------- 

106 cluster_space 

107 Cluster space implicitly defining allowed species. 

108 structure 

109 Specific supercell (consistent with cluster_space) whose 

110 allowed species are to be determined. 

111 """ 

112 

113 primitive_structure = Structure.from_atoms(cluster_space.primitive_structure) 

114 chemical_symbols = cluster_space.chemical_symbols 

115 

116 if len(chemical_symbols) == 1: 

117 # If the allowed species are the same for all sites no loop is 

118 # required 

119 allowed_species = {site: chemical_symbols[0] for 

120 site in self._sites.keys()} 

121 else: 

122 # Loop over the lattice sites to find the allowed species 

123 allowed_species = {} 

124 for site, indices in self._sites.items(): 

125 allowed_species[site] = None 

126 positions = structure.positions[np.array(indices)] 

127 lattice_sites = primitive_structure.find_lattice_sites_by_positions( 

128 positions=positions, 

129 fractional_position_tolerance=cluster_space.fractional_position_tolerance) 

130 for k, lattice_site in enumerate(lattice_sites): 

131 species = chemical_symbols[lattice_site.index] 

132 # check that the allowed species are equal for all sites 

133 if allowed_species[site] is not None and \ 

134 species != allowed_species[site]: 

135 raise Exception('The allowed species {} for the site' 

136 ' with index {} differs from the' 

137 ' result {} for the previous index' 

138 ' ({})!'.format(species, indices[k], 

139 allowed_species[site], 

140 indices[k-1])) 

141 allowed_species[site] = species 

142 

143 self._allowed_species = allowed_species 

144 

145 def get_observable(self, structure: Atoms) -> Dict[str, List[float]]: 

146 """ 

147 Returns the site occupation factors for a given atomic configuration. 

148 

149 Parameters 

150 ---------- 

151 structure 

152 Input atomic structure. 

153 """ 

154 

155 chemical_symbols = np.array(structure.get_chemical_symbols()) 

156 sofs = {} 

157 for site, indices in self._sites.items(): 

158 counts = {species: 0 for species in self._allowed_species[site]} 

159 symbols, sym_counts = np.unique(chemical_symbols[indices], 

160 return_counts=True) 

161 for sym, count in zip(symbols, sym_counts): 

162 counts[sym] += count 

163 

164 for species in counts.keys(): 

165 key = 'sof_{}_{}'.format(site, species) 

166 sofs[key] = float(counts[species]) / len(indices) 

167 

168 return sofs