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

1from math import inf 

2import numpy as np 

3from typing import List, Dict 

4 

5from ase import Atoms 

6from ase.data import chemical_symbols as periodic_table 

7from .. import ClusterExpansion 

8from ..core.local_orbit_list_generator import LocalOrbitListGenerator 

9from ..core.structure import Structure 

10from .variable_transformation import transform_parameters 

11from ..input_output.logging_tools import logger 

12from pkg_resources import VersionConflict 

13 

14try: 

15 import mip 

16 from mip.constants import BINARY, INTEGER 

17 from distutils.version import LooseVersion 

18 

19 if LooseVersion(mip.constants.VERSION) < '1.6.3': 

20 raise VersionConflict('Python-MIP version 1.6.3 or later is required in order to use the ' 

21 'ground state finder.') 

22except ImportError: 

23 raise ImportError('Python-MIP (https://python-mip.readthedocs.io/en/latest/) is required in ' 

24 'order to use the ground state finder.') 

25 

26 

27class GroundStateFinder: 

28 """ 

29 This class provides functionality for determining the ground states 

30 using a binary cluster expansion. This is efficiently achieved through the 

31 use of mixed integer programming (MIP) as developed by Larsen *et al.* in 

32 `Phys. Rev. Lett. 120, 256101 (2018) 

33 <https://doi.org/10.1103/PhysRevLett.120.256101>`_. 

34 

35 This class relies on the `Python-MIP package 

36 <https://python-mip.readthedocs.io>`_. Python-MIP can be used together 

37 with `Gurobi <https://www.gurobi.com/>`_, which is not open source 

38 but issues academic licenses free of charge. Pleaase note that 

39 Gurobi needs to be installed separately. The `GroundStateFinder` works 

40 also without Gurobi, but if performance is critical, Gurobi is highly 

41 recommended. 

42 

43 Warning 

44 ------- 

45 In order to be able to use Gurobi with python-mip one must ensure that 

46 `GUROBI_HOME` should point to the installation directory 

47 (``<installdir>``):: 

48 

49 export GUROBI_HOME=<installdir> 

50 

51 Note 

52 ---- 

53 The current implementation only works for binary systems. 

54 

55 

56 Parameters 

57 ---------- 

58 cluster_expansion : ClusterExpansion 

59 cluster expansion for which to find ground states 

60 structure : Atoms 

61 atomic configuration 

62 solver_name : str, optional 

63 'gurobi', alternatively 'grb', or 'cbc', searches for available 

64 solvers if not informed 

65 verbose : bool, optional 

66 whether to display solver messages on the screen 

67 (default: True) 

68 

69 

70 Example 

71 ------- 

72 The following snippet illustrates how to determine the ground state for a 

73 Au-Ag alloy. Here, the parameters of the cluster 

74 expansion are set to emulate a simple Ising model in order to obtain an 

75 example that can be run without modification. In practice, one should of 

76 course use a proper cluster expansion:: 

77 

78 >>> from ase.build import bulk 

79 >>> from icet import ClusterExpansion, ClusterSpace 

80 

81 >>> # prepare cluster expansion 

82 >>> # the setup emulates a second nearest-neighbor (NN) Ising model 

83 >>> # (zerolet and singlet parameters are zero; only first and second neighbor 

84 >>> # pairs are included) 

85 >>> prim = bulk('Au') 

86 >>> chemical_symbols = ['Ag', 'Au'] 

87 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=chemical_symbols) 

88 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02]) 

89 

90 >>> # prepare initial configuration 

91 >>> structure = prim.repeat(3) 

92 

93 >>> # set up the ground state finder and calculate the ground state energy 

94 >>> gsf = GroundStateFinder(ce, structure) 

95 >>> ground_state = gsf.get_ground_state({'Ag': 5}) 

96 >>> print('Ground state energy:', ce.predict(ground_state)) 

97 """ 

98 

99 def __init__(self, 

100 cluster_expansion: ClusterExpansion, 

101 structure: Atoms, 

102 solver_name: str = None, 

103 verbose: bool = True) -> None: 

104 # Check that there is only one active sublattice 

105 self._cluster_expansion = cluster_expansion 

106 self._fractional_position_tolerance = cluster_expansion.fractional_position_tolerance 

107 self.structure = structure 

108 cluster_space = self._cluster_expansion.get_cluster_space_copy() 

109 primitive_structure = cluster_space.primitive_structure 

110 self._active_sublattices = cluster_space.get_sublattices(structure).active_sublattices 

111 

112 # Check that there are no more than two allowed species 

113 active_species = [set(subl.chemical_symbols) for subl in self._active_sublattices] 

114 if any(len(species) > 2 for species in active_species): 

115 raise NotImplementedError('Currently, systems with more than two allowed species on ' 

116 'any sublattice are not supported.') 

117 self._active_species = active_species 

118 

119 # Define cluster functions for elements 

120 self._reverse_id_maps = [] 

121 for species in active_species: 

122 for species_map in cluster_space.species_maps: 122 ↛ 121line 122 didn't jump to line 121, because the loop on line 122 didn't complete

123 symbols = [periodic_table[n] for n in species_map] 

124 if set(symbols) == species: 

125 reverse_id_map = {1 - species_map[n]: periodic_table[n] for n in species_map} 

126 self._reverse_id_maps.append(reverse_id_map) 

127 break 

128 self._count_symbols = [reverse_id_map[1] for reverse_id_map in self._reverse_id_maps] 

129 

130 # Generate full orbit list 

131 self._orbit_list = cluster_space.orbit_list 

132 lolg = LocalOrbitListGenerator( 

133 orbit_list=self._orbit_list, 

134 structure=Structure.from_atoms(primitive_structure), 

135 fractional_position_tolerance=self._fractional_position_tolerance) 

136 self._full_orbit_list = lolg.generate_full_orbit_list() 

137 

138 # Transform the parameters 

139 binary_parameters = transform_parameters(primitive_structure, 

140 self._full_orbit_list, 

141 self._cluster_expansion.parameters) 

142 self._transformed_parameters = binary_parameters 

143 

144 # Build model 

145 if solver_name is None: 

146 solver_name = '' 

147 self._model = self._build_model(structure, solver_name, verbose) 

148 

149 # Properties that are defined when searching for a ground state 

150 self._optimization_status = None 

151 

152 def _build_model(self, 

153 structure: Atoms, 

154 solver_name: str, 

155 verbose: bool) -> mip.Model: 

156 """ 

157 Build a Python-MIP model based on the provided structure 

158 

159 Parameters 

160 ---------- 

161 structure 

162 atomic configuration 

163 solver_name 

164 'gurobi', alternatively 'grb', or 'cbc', searches for 

165 available solvers if not informed 

166 verbose 

167 whether to display solver messages on the screen 

168 """ 

169 

170 # Create cluster maps 

171 self._create_cluster_maps(structure) 

172 

173 # Initiate MIP model 

174 model = mip.Model('CE', solver_name=solver_name) 

175 model.solver.set_mip_gap(0) # avoid stopping prematurely 

176 model.solver.set_emphasis(2) # focus on finding optimal solution 

177 model.preprocess = 2 # maximum preprocessing 

178 

179 # Set verbosity 

180 model.verbose = int(verbose) 

181 

182 # Spin variables (remapped) for all atoms in the structure 

183 xs = {i: model.add_var(name='atom_{}'.format(i), var_type=BINARY) 

184 for subl in self._active_sublattices for i in subl.indices} 

185 ys = [model.add_var(name='cluster_{}'.format(i), var_type=BINARY) 

186 for i in range(len(self._cluster_to_orbit_map))] 

187 

188 # The objective function is added to 'model' first 

189 model.objective = mip.minimize(mip.xsum(self._get_total_energy(ys))) 

190 

191 # Connect cluster variables to spin variables with cluster constraints 

192 # TODO: don't create cluster constraints for singlets 

193 constraint_count = 0 

194 for i, cluster in enumerate(self._cluster_to_sites_map): 

195 orbit = self._cluster_to_orbit_map[i] 

196 parameter = self._transformed_parameters[orbit + 1] 

197 assert parameter != 0 

198 

199 if len(cluster) < 2 or parameter < 0: # no "downwards" pressure 

200 for atom in cluster: 

201 model.add_constr(ys[i] <= xs[atom], 

202 'Decoration -> cluster {}'.format(constraint_count)) 

203 constraint_count = constraint_count + 1 

204 

205 if len(cluster) < 2 or parameter > 0: # no "upwards" pressure 

206 model.add_constr(ys[i] >= 1 - len(cluster) + 

207 mip.xsum(xs[atom] 

208 for atom in cluster), 

209 'Decoration -> cluster {}'.format(constraint_count)) 

210 constraint_count = constraint_count + 1 

211 

212 for sym, subl in zip(self._count_symbols, self._active_sublattices): 

213 # Create slack variable 

214 slack = model.add_var(name='slackvar_{}'.format(sym), var_type=INTEGER, 

215 lb=0, ub=len(subl.indices)) 

216 

217 # Add slack constraint 

218 model.add_constr(slack <= -1, name='{} slack'.format(sym)) 

219 

220 # Set species constraint 

221 model.add_constr(mip.xsum([xs[i] for i in subl.indices]) + slack == -1, 

222 name='{} count'.format(sym)) 

223 

224 # Update the model so that variables and constraints can be queried 

225 if model.solver_name.upper() in ['GRB', 'GUROBI']: 225 ↛ 226line 225 didn't jump to line 226, because the condition on line 225 was never true

226 model.solver.update() 

227 return model 

228 

229 def _create_cluster_maps(self, structure: Atoms) -> None: 

230 """ 

231 Create maps that include information regarding which sites and orbits 

232 are associated with each cluster as well as the number of clusters per 

233 orbit 

234 

235 Parameters 

236 ---------- 

237 structure 

238 atomic configuration 

239 """ 

240 # Generate full orbit list 

241 lolg = LocalOrbitListGenerator( 

242 orbit_list=self._orbit_list, 

243 structure=Structure.from_atoms(structure), 

244 fractional_position_tolerance=self._fractional_position_tolerance) 

245 full_orbit_list = lolg.generate_full_orbit_list() 

246 

247 # Create maps of site indices and orbits for all clusters 

248 cluster_to_sites_map = [] 

249 cluster_to_orbit_map = [] 

250 for orb_index in range(len(full_orbit_list)): 

251 

252 equivalent_clusters = full_orbit_list.get_orbit(orb_index).equivalent_clusters 

253 

254 # Determine the sites and the orbit associated with each cluster 

255 for cluster in equivalent_clusters: 

256 

257 # Do not include clusters for which the parameter is 0 

258 parameter = self._transformed_parameters[orb_index + 1] 

259 if parameter == 0: 

260 continue 

261 

262 # Add the the list of sites and the orbit to the respective cluster maps 

263 cluster_sites = [site.index for site in cluster] 

264 cluster_to_sites_map.append(cluster_sites) 

265 cluster_to_orbit_map.append(orb_index) 

266 

267 # calculate the number of clusters per orbit 

268 nclusters_per_orbit = [cluster_to_orbit_map.count(i) for i in 

269 range(cluster_to_orbit_map[-1] + 1)] 

270 nclusters_per_orbit = [1] + nclusters_per_orbit 

271 

272 self._cluster_to_sites_map = cluster_to_sites_map 

273 self._cluster_to_orbit_map = cluster_to_orbit_map 

274 self._nclusters_per_orbit = nclusters_per_orbit 

275 

276 def _get_total_energy(self, cluster_instance_activities: List[int]) -> List[float]: 

277 """ 

278 Calculates the total energy using the expression based on binary 

279 variables 

280 

281 .. math:: 

282 

283 H({\\boldsymbol x}, {\\boldsymbol E})=E_0+ 

284 \\sum\\limits_j\\sum\\limits_{{\\boldsymbol c} 

285 \\in{\\boldsymbol C}_j}E_jy_{{\\boldsymbol c}}, 

286 

287 where (:math:`y_{{\\boldsymbol c}}= 

288 \\prod\\limits_{i\\in{\\boldsymbol c}}x_i`). 

289 

290 Parameters 

291 ---------- 

292 cluster_instance_activities 

293 list of cluster instance activities, (:math:`y_{{\\boldsymbol c}}`) 

294 """ 

295 

296 E = [0.0 for _ in self._transformed_parameters] 

297 for i in range(len(cluster_instance_activities)): 

298 orbit = self._cluster_to_orbit_map[i] 

299 E[orbit + 1] = E[orbit + 1] + cluster_instance_activities[i] 

300 E[0] = 1 

301 

302 E = [0.0 if np.isclose(self._transformed_parameters[orbit], 0.0) else 

303 E[orbit] * self._transformed_parameters[orbit] / self._nclusters_per_orbit[orbit] 

304 for orbit in range(len(self._transformed_parameters))] 

305 return E 

306 

307 def get_ground_state(self, 

308 species_count: Dict[str, int] = None, 

309 max_seconds: float = inf, 

310 threads: int = 0) -> Atoms: 

311 """ 

312 Finds the ground state for a given structure and species count, which 

313 refers to the `count_species`, if provided when initializing the 

314 instance of this class, or the first species in the list of chemical 

315 symbols for the active sublattice. 

316 

317 Parameters 

318 ---------- 

319 species_count 

320 dictionary with count for one of the species on each active 

321 sublattice. If no count is provided for a sublattice, the 

322 concentration is allowed to vary. 

323 max_seconds 

324 maximum runtime in seconds (default: inf) 

325 threads 

326 number of threads to be used when solving the problem, given that a 

327 positive integer has been provided. If set to 0 the solver default 

328 configuration is used while -1 corresponds to all available 

329 processing cores. 

330 """ 

331 if species_count is None: 

332 species_count = {} 

333 

334 # Check that the species_count is consistent with the cluster space 

335 all_active_species = set.union(*self._active_species) 

336 for symbol in species_count: 

337 if symbol not in all_active_species: 

338 raise ValueError('The species {} is not present on any of the active sublattices' 

339 ' ({})'.format(symbol, self._active_species)) 

340 

341 # The model is solved using python-MIPs choice of solver, which is 

342 # Gurobi, if available, and COIN-OR Branch-and-Cut, otherwise. 

343 model = self._model 

344 

345 # Update the species counts 

346 for i, species in enumerate(self._active_species): 

347 count_symbol = self._count_symbols[i] 

348 max_count = len(self._active_sublattices[i].indices) 

349 

350 symbols_to_add = set.intersection(set(species_count), set(species)) 

351 if len(symbols_to_add) > 1: 

352 raise ValueError('Provide counts for at most one of the species on each active ' 

353 'sublattice ({}), not {}!'.format(self._active_species, 

354 list(species_count))) 

355 elif len(symbols_to_add) == 1: 

356 sym = symbols_to_add.pop() 

357 count = species_count[sym] 

358 if count < 0 or count > max_count: 

359 raise ValueError('The count for species {} ({}) must be a positive integer and' 

360 ' cannot exceed the number of sites on the active sublattice ' 

361 '({})'.format(sym, count, max_count)) 

362 if sym == count_symbol: 

363 xcount = count 

364 else: 

365 xcount = max_count - count 

366 

367 max_slack = 0 

368 else: 

369 xcount = max_slack = max_count 

370 

371 model.constr_by_name('{} count'.format(count_symbol)).rhs = xcount 

372 model.constr_by_name('{} slack'.format(count_symbol)).rhs = max_slack 

373 

374 # Set the number of threads 

375 model.threads = threads 

376 

377 # Optimize the model 

378 self._optimization_status = model.optimize(max_seconds=max_seconds) 

379 

380 # The status of the solution is printed to the screen 

381 if str(self._optimization_status) != 'OptimizationStatus.OPTIMAL': 

382 if str(self._optimization_status) == 'OptimizationStatus.FEASIBLE': 

383 logger.warning('Solution optimality not proven.') 

384 else: 

385 raise Exception('Optimization failed ({0})'.format(str(self._optimization_status))) 

386 

387 # Each of the variables is printed with it's resolved optimum value 

388 gs = self.structure.copy() 

389 

390 active_index_to_sublattice_map = {i: j for j, subl in enumerate(self._active_sublattices) 

391 for i in subl.indices} 

392 for v in model.vars: 

393 if 'atom' in v.name: 

394 index = int(v.name.split('_')[-1]) 

395 sublattice_index = active_index_to_sublattice_map[index] 

396 gs[index].symbol = self._reverse_id_maps[sublattice_index][int(v.x)] 

397 

398 # Assert that the solution agrees with the prediction 

399 prediction = self._cluster_expansion.predict(gs) 

400 assert abs(model.objective_value - prediction) < 1e-6 

401 return gs 

402 

403 @property 

404 def optimization_status(self) -> mip.OptimizationStatus: 

405 """Optimization status""" 

406 return self._optimization_status 

407 

408 @property 

409 def model(self) -> mip.Model: 

410 """Python-MIP model""" 

411 return self._model.copy()