Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

from ase import Atoms 

from icet import ClusterSpace 

from icet.core.structure import Structure 

from mchammer.observers.base_observer import BaseObserver 

from typing import List, Dict 

import numpy as np 

 

 

class SiteOccupancyObserver(BaseObserver): 

""" 

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

 

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

trajectory sampled by a Monte Carlo (MC) simulation. 

 

Parameters 

---------- 

cluster_space : icet.ClusterSpace 

cluster space from which the allowed species are extracted 

 

sites : dict(str, list(int)) 

dictionary containing lists of sites that are to be considered, 

which keys will be taken as the names of the sites 

 

supercell : ase.Atoms 

a typical supercell, which is used to determine the allowed species 

 

interval : int 

observation interval during the Monte Carlo simulation 

 

Attributes 

---------- 

tag : str 

name of observer 

 

interval : int 

observation interval 

 

Example 

------- 

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

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

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

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

cluster expansion is used with slightly modified surface interactions in 

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

one should of course use a proper cluster expansion:: 

 

from ase.build import fcc111 

from icet import ClusterExpansion, ClusterSpace 

from mchammer.calculators import ClusterExpansionCalculator 

from mchammer.ensembles import CanonicalEnsemble 

from mchammer.observers import SiteOccupancyObserver 

 

# prepare reference structure 

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

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

prim.wrap() 

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

 

# prepare cluster expansion 

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

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

params[1] = 0.01 

params[6] = 0.12 

ce = ClusterExpansion(cs, params) 

print(ce) 

 

# prepare initial configuration based on a 2x2 supercell 

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

for k in range(20): 

structure[k].symbol = 'Ag' 

 

# set up MC simulation 

calc = ClusterExpansionCalculator(structure, ce) 

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

data_container='myrun_sof.dc') 

 

# set up observer and attach it to the MC simulation 

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

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

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

mc.attach_observer(sof) 

 

# run 1000 trial steps 

mc.run(1000) 

 

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

container:: 

 

print(mc.data_container.data) 

""" 

 

def __init__(self, cluster_space: ClusterSpace, 

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

supercell: Atoms, 

interval: int = None) -> None: 

super().__init__(interval=interval, return_type=dict, 

tag='SiteOccupancyObserver') 

 

self._sites = {site: sorted(indices) 

for site, indices in sites.items()} 

 

self._set_allowed_species(cluster_space, supercell) 

 

def _set_allowed_species(self, 

cluster_space: ClusterSpace, 

supercell: Atoms): 

""" 

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

 

Parameters 

---------- 

cluster_space 

Cluster space implicitly defining allowed species 

supercell 

Specific supercell (consistent with cluster_space) whose 

allowed species are to be determined 

""" 

 

primitive_structure = Structure.from_atoms( 

cluster_space.primitive_structure) 

chemical_symbols = cluster_space.get_chemical_symbols() 

 

if len(chemical_symbols) == 1: 

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

# required 

allowed_species = {site: chemical_symbols[0] for 

site in self._sites.keys()} 

else: 

# Loop over the lattice sites to find the allowed species 

allowed_species = {} 

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

allowed_species[site] = None 

positions = supercell.get_positions()[np.array(indices)] 

lattice_sites =\ 

primitive_structure.find_lattice_sites_by_positions( 

positions) 

for l, lattice_site in enumerate(lattice_sites): 

species = chemical_symbols[lattice_site.index] 

# check that the allowed species are equal for all sites 

if allowed_species[site] is not None and\ 

species != allowed_species[site]: 

raise Exception("The allowed species {} for the site" 

" with index {} differs from the" 

" result {} for the previous index" 

" ({})!".format(species, indices[l], 

allowed_species[site], 

indices[l-1])) 

allowed_species[site] = species 

 

self._allowed_species = allowed_species 

 

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

""" 

Returns the site occupation factors for a given atomic configuration. 

 

Parameters 

---------- 

structure 

input atomic structure. 

""" 

 

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

sofs = {} 

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

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

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

return_counts=True) 

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

counts[sym] += count 

 

for species in counts.keys(): 

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

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

 

return sofs