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

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

from abc import abstractproperty 

from typing import Dict, List 

 

import numpy as np 

 

from ase import Atoms 

from ase.units import kB 

from ase.data import chemical_symbols 

 

from .. import DataContainer 

from ..calculators.base_calculator import BaseCalculator 

from .base_ensemble import BaseEnsemble 

 

 

class ThermodynamicBaseEnsemble(BaseEnsemble): 

""" 

Thermodynamic base ensemble class. 

 

Parameters 

---------- 

structure : :class:`Atoms <ase.Atoms>` 

atomic configuration to be used in the Monte Carlo simulation; 

also defines the initial occupation vector 

calculator : :class:`BaseCalculator <mchammer.calculators.ClusterExpansionCalculator>` 

calculator to be used for calculating the potential changes 

that enter the evaluation of the Metropolis criterion 

boltzmann_constant : float 

Boltzmann constant :math:`k_B` in appropriate 

units, i.e. units that are consistent 

with the underlying cluster expansion 

and the temperature units [default: eV/K] 

user_tag : str 

human-readable tag for ensemble [default: None] 

data_container : str 

name of file the data container associated with the ensemble 

will be written to; if the file exists it will be read, the 

data container will be appended, and the file will be 

updated/overwritten 

random_seed : int 

seed for the random number generator used in the Monte Carlo 

simulation 

ensemble_data_write_interval : int 

interval at which data is written to the data container; this 

includes for example the current value of the calculator 

(i.e. usually the energy) as well as ensembles specific fields 

such as temperature or the number of structure of different species 

data_container_write_period : float 

period in units of seconds at which the data container is 

written to file; writing periodically to file provides both 

a way to examine the progress of the simulation and to back up 

the data [default: np.inf] 

trajectory_write_interval : int 

interval at which the current occupation vector of the atomic 

configuration is written to the data container. 

""" 

 

def __init__(self, 

structure: Atoms, 

calculator: BaseCalculator, 

user_tag: str = None, 

boltzmann_constant: float = kB, 

data_container: DataContainer = None, 

random_seed: int = None, 

data_container_write_period: float = np.inf, 

ensemble_data_write_interval: int = None, 

trajectory_write_interval: int = None) -> None: 

 

self._boltzmann_constant = boltzmann_constant 

 

super().__init__( 

structure=structure, calculator=calculator, user_tag=user_tag, 

data_container=data_container, 

random_seed=random_seed, 

data_container_write_period=data_container_write_period, 

ensemble_data_write_interval=ensemble_data_write_interval, 

trajectory_write_interval=trajectory_write_interval) 

 

@abstractproperty 

@property 

def temperature(self) -> float: 

pass 

 

@property 

def boltzmann_constant(self) -> float: 

""" Boltzmann constant :math:`k_B` (see parameters section above) """ 

return self._boltzmann_constant 

 

def _acceptance_condition(self, potential_diff: float) -> bool: 

""" 

Evaluates Metropolis acceptance criterion. 

 

Parameters 

---------- 

potential_diff 

change in the thermodynamic potential associated 

with the trial step 

""" 

if potential_diff <= 0: 

return True 

elif self.temperature <= 1e-16: 

return False 

else: 

p = np.exp(-potential_diff / (self.boltzmann_constant * self.temperature)) 

return p > self._next_random_number() 

 

def do_canonical_swap(self, sublattice_index: int, allowed_species: List[int] = None) -> int: 

""" Carries out one Monte Carlo trial step. 

 

Parameters 

--------- 

sublattice_index 

the sublattice the swap will act on 

allowed_species 

list of atomic numbers for allowed species 

 

Returns 

------- 

Returns 1 or 0 depending on if trial move was accepted or rejected 

""" 

sites, species = self.configuration.get_swapped_state(sublattice_index, allowed_species) 

potential_diff = self._get_property_change(sites, species) 

 

if self._acceptance_condition(potential_diff): 

self.update_occupations(sites, species) 

return 1 

return 0 

 

def do_sgc_flip(self, chemical_potentials: Dict[int, float], sublattice_index: int, 

allowed_species: List[int] = None) -> int: 

""" Carries out one Monte Carlo trial step. 

 

Parameters 

--------- 

chemical_potentials 

chemical potentials used to calculate the potential 

difference 

sublattice_index 

the sublattice the flip will act on 

allowed_species 

list of atomic numbers for allowed species 

 

Returns 

------- 

Returns 1 or 0 depending on if trial move was accepted or rejected 

""" 

index, species = self.configuration.get_flip_state(sublattice_index, allowed_species) 

potential_diff = self._get_property_change([index], [species]) 

 

# change in chemical potential 

old_species = self.configuration.occupations[index] 

chemical_potential_diff = chemical_potentials[old_species] - chemical_potentials[species] 

potential_diff += chemical_potential_diff 

 

if self._acceptance_condition(potential_diff): 

self.update_occupations([index], [species]) 

return 1 

return 0 

 

def do_vcsgc_flip(self, phis: Dict[int, float], kappa: float, sublattice_index: int, 

allowed_species: List[int] = None) -> int: 

"""Carries out one Monte Carlo trial step. 

 

Parameters 

---------- 

phis 

average constraint parameters 

kappa 

parameter that constrains the variance of the concentration 

sublattice_index 

the sublattice the flip will act on 

allowed_species 

list of atomic numbers for allowed species 

 

Returns 

------- 

Returns 1 or 0 depending on if trial move was accepted or rejected 

""" 

index, new_species = self.configuration.get_flip_state( 

sublattice_index, allowed_species) 

old_species = self.configuration.occupations[index] 

 

# Calculate difference in VCSGC thermodynamic potential. 

# Note that this assumes that only one atom was flipped. 

sl_occupations = self.configuration.get_occupations_on_sublattice(sublattice_index) 

N = len(sl_occupations) 

potential_diff = 1.0 # dN 

for species in phis: 

if species == old_species: 

factor = -1 

190 ↛ 193line 190 didn't jump to line 193, because the condition on line 190 was never false elif species == new_species: 

factor = 1 

else: 

continue 

potential_diff += factor * (N * phis[species] + 2 * sl_occupations.count(species)) 

potential_diff *= kappa * self.boltzmann_constant * self.temperature / N 

potential_diff += self._get_property_change([index], [new_species]) 

 

if self._acceptance_condition(potential_diff): 

self.update_occupations([index], [new_species]) 

return 1 

return 0 

 

def _get_swap_sublattice_probabilities(self) -> List[float]: 

""" Returns sublattice probabilities suitable for swaps.""" 

sublattice_probabilities = [] 

for i, sl in enumerate(self.sublattices): 

if self.configuration.is_swap_possible(i): 

sublattice_probabilities.append(len(sl.indices)) 

else: 

sublattice_probabilities.append(0) 

norm = sum(sublattice_probabilities) 

if norm == 0: 

raise ValueError('No canonical swaps are possible on any of the active sublattices.') 

sublattice_probabilities = [p / norm for p in sublattice_probabilities] 

return sublattice_probabilities 

 

def _get_flip_sublattice_probabilities(self) -> List[float]: 

"""Returns the default sublattice probability which is based on 

the sizes of a sublattice. 

""" 

probability_distribution = [] 

for i, sl in enumerate(self.sublattices): 

if len(sl.chemical_symbols) > 1: 

probability_distribution.append(len(sl.indices)) 

else: 

probability_distribution.append(0) 

 

norm = sum(probability_distribution) 

probability_distribution = [p / norm for p in probability_distribution] 

return probability_distribution 

 

def _get_vcsgc_free_energy_derivatives(self, phis: Dict[int, str], kappa: float, 

sublattice_index: int = None) -> Dict: 

""" 

Returns a dict with the free energy derivatives. 

 

Parameters 

---------- 

phis 

average constraint parameters 

kappa 

parameter that constrains the variance of the concentration 

sublattice_index 

sublattice index 

""" 

data = {} 

 

for atnum in phis: 

for i, sublattice in enumerate(self.sublattices): 

if sublattice_index is not None and i != sublattice_index: 

continue 

if len(sublattice.chemical_symbols) > 0 and atnum in sublattice.atomic_numbers: 

N = len(sublattice.indices) 

sl_occupations = self.configuration.get_occupations_on_sublattice(i) 

concentration = sl_occupations.count(atnum) / N 

data['free_energy_derivative_{}'.format(chemical_symbols[atnum])] \ 

= kappa * self.boltzmann_constant * self.temperature * \ 

(- 2 * concentration - phis[atnum]) 

 

return data 

 

def _get_species_counts(self) -> Dict: 

""" 

Returns a dict with the species counts. 

""" 

data = {} 

structure = self.configuration.structure 

unique, counts = np.unique(structure.numbers, return_counts=True) 

for sl in self.sublattices: 

for symbol in sl.chemical_symbols: 

data['{}_count'.format(symbol)] = 0 

for atnum, count in zip(unique, counts): 

data['{}_count'.format(chemical_symbols[atnum])] = count 

 

return data