from math import inf
import numpy as np
from typing import List, Dict
from ase import Atoms
from ase.data import chemical_symbols as periodic_table
from .. import ClusterExpansion
from ..core.local_orbit_list_generator import LocalOrbitListGenerator
from ..core.structure import Structure
from .variable_transformation import transform_parameters
from ..input_output.logging_tools import logger
try:
import mip
from mip.constants import BINARY, INTEGER
except ImportError:
raise ImportError('Python-MIP (https://python-mip.readthedocs.io/en/latest/) is required in '
'order to use the ground state finder.')
[docs]
class GroundStateFinder:
"""
This class provides functionality for determining the ground states
using a binary cluster expansion. This is efficiently achieved through the
use of mixed integer programming (MIP) as developed by Larsen *et al.* in
`Phys. Rev. Lett. 120, 256101 (2018)
<https://doi.org/10.1103/PhysRevLett.120.256101>`_.
This class relies on the `Python-MIP package
<https://python-mip.readthedocs.io>`_. Python-MIP can be used together
with `Gurobi <https://www.gurobi.com/>`_, which is not open source
but issues academic licenses free of charge. Pleaase note that
Gurobi needs to be installed separately. The :class:`GroundStateFinder` works
also without Gurobi, but if performance is critical, Gurobi is highly
recommended.
Warning
-------
In order to be able to use Gurobi with python-mip one must ensure that
`GUROBI_HOME` should point to the installation directory
(``<installdir>``)::
export GUROBI_HOME=<installdir>
Note
----
The current implementation only works for binary systems.
Parameters
----------
cluster_expansion
Cluster expansion for which to find ground states.
structure
Atomic configuration.
solver_name
``'gurobi'``, ``'grb'`` or ``'cbc'``. Searches for available
solvers if no value is provided.
verbose
If ``True`` print solver messages to stdout.
Example
-------
The following snippet illustrates how to determine the ground state for a
Au-Ag alloy. Here, the parameters of the cluster
expansion are set to emulate a simple Ising model in order to obtain an
example that can be run without modification. In practice, one should of
course use a proper cluster expansion::
>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet parameters are zero; only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> chemical_symbols = ['Ag', 'Au']
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=chemical_symbols)
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
>>> # prepare initial configuration
>>> structure = prim.repeat(3)
>>> # set up the ground state finder and calculate the ground state energy
>>> gsf = GroundStateFinder(ce, structure)
>>> ground_state = gsf.get_ground_state({'Ag': 5})
>>> print('Ground state energy:', ce.predict(ground_state))
"""
def __init__(self,
cluster_expansion: ClusterExpansion,
structure: Atoms,
solver_name: str = None,
verbose: bool = True) -> None:
# Check that there is only one active sublattice
self._cluster_expansion = cluster_expansion
self._fractional_position_tolerance = cluster_expansion.fractional_position_tolerance
self.structure = structure
cluster_space = self._cluster_expansion.get_cluster_space_copy()
primitive_structure = cluster_space.primitive_structure
self._active_sublattices = cluster_space.get_sublattices(structure).active_sublattices
# Check that there are no more than two allowed species
active_species = [set(subl.chemical_symbols) for subl in self._active_sublattices]
if any(len(species) > 2 for species in active_species):
raise NotImplementedError('Currently, systems with more than two allowed species on '
'any sublattice are not supported.')
# Check that there are no merged orbits
if any(['merge' in elem for elem in cluster_space._pruning_history]):
raise NotImplementedError('Currently, systems with merged orbits are not supported.')
self._active_species = active_species
# Define cluster functions for elements
self._reverse_id_maps = []
for species in active_species:
for species_map in cluster_space.species_maps:
symbols = [periodic_table[n] for n in species_map]
if set(symbols) == species:
reverse_id_map = {1 - species_map[n]: periodic_table[n] for n in species_map}
self._reverse_id_maps.append(reverse_id_map)
break
self._count_symbols = [reverse_id_map[1] for reverse_id_map in self._reverse_id_maps]
# Generate full orbit list
self._orbit_list = cluster_space.orbit_list
lolg = LocalOrbitListGenerator(
orbit_list=self._orbit_list,
structure=Structure.from_atoms(primitive_structure),
fractional_position_tolerance=self._fractional_position_tolerance)
self._full_orbit_list = lolg.generate_full_orbit_list()
# Transform the parameters
binary_parameters = transform_parameters(primitive_structure,
self._full_orbit_list,
self._cluster_expansion.parameters)
self._transformed_parameters = binary_parameters
# Build model
if solver_name is None:
solver_name = ''
self._model = self._build_model(structure, solver_name, verbose)
# Properties that are defined when searching for a ground state
self._optimization_status = None
def _build_model(self,
structure: Atoms,
solver_name: str,
verbose: bool) -> mip.Model:
"""
Build a Python-MIP model based on the provided structure.
Parameters
----------
structure
Atomic configuration.
solver_name
``'gurobi'``, ``'grb'`` or ``'cbc'``. Searches for available
solvers if no value is provided.
verbose
If ``True`` print solver messages to stdout.
"""
# Create cluster maps
self._create_cluster_maps(structure)
# Initiate MIP model
model = mip.Model('CE', solver_name=solver_name)
model.solver.set_mip_gap(0) # avoid stopping prematurely
model.solver.set_emphasis(2) # focus on finding optimal solution
model.preprocess = 2 # maximum preprocessing
# Set verbosity
model.verbose = int(verbose)
# Spin variables (remapped) for all atoms in the structure
xs = {i: model.add_var(name='atom_{}'.format(i), var_type=BINARY)
for subl in self._active_sublattices for i in subl.indices}
ys = [model.add_var(name='cluster_{}'.format(i), var_type=BINARY)
for i in range(len(self._cluster_to_orbit_map))]
# The objective function is added to 'model' first
model.objective = mip.minimize(mip.xsum(self._get_total_energy(ys)))
# Connect cluster variables to spin variables with cluster constraints
# TODO: don't create cluster constraints for singlets
constraint_count = 0
for i, cluster in enumerate(self._cluster_to_sites_map):
orbit = self._cluster_to_orbit_map[i]
parameter = self._transformed_parameters[orbit + 1]
assert parameter != 0
if len(cluster) < 2 or parameter < 0: # no "downwards" pressure
for atom in cluster:
model.add_constr(ys[i] <= xs[atom],
'Decoration -> cluster {}'.format(constraint_count))
constraint_count = constraint_count + 1
if len(cluster) < 2 or parameter > 0: # no "upwards" pressure
model.add_constr(ys[i] >= 1 - len(cluster) +
mip.xsum(xs[atom]
for atom in cluster),
'Decoration -> cluster {}'.format(constraint_count))
constraint_count = constraint_count + 1
for sym, subl in zip(self._count_symbols, self._active_sublattices):
# Create slack variable
slack = model.add_var(name='slackvar_{}'.format(sym), var_type=INTEGER,
lb=0, ub=len(subl.indices))
# Add slack constraint
model.add_constr(slack <= -1, name='{} slack'.format(sym))
# Set species constraint
model.add_constr(mip.xsum([xs[i] for i in subl.indices]) + slack == -1,
name='{} count'.format(sym))
# Update the model so that variables and constraints can be queried
if model.solver_name.upper() in ['GRB', 'GUROBI']:
model.solver.update()
return model
def _create_cluster_maps(self, structure: Atoms) -> None:
"""
Create maps that include information regarding which sites and orbits
are associated with each cluster as well as the number of clusters per
orbit.
Parameters
----------
structure
Atomic configuration.
"""
# Generate full orbit list
lolg = LocalOrbitListGenerator(
orbit_list=self._orbit_list,
structure=Structure.from_atoms(structure),
fractional_position_tolerance=self._fractional_position_tolerance)
full_orbit_list = lolg.generate_full_orbit_list()
# Create maps of site indices and orbits for all clusters
cluster_to_sites_map = []
cluster_to_orbit_map = []
for orb_index in range(len(full_orbit_list)):
clusters = full_orbit_list.get_orbit(orb_index).clusters
# Determine the sites and the orbit associated with each cluster
for cluster in clusters:
# Do not include clusters for which the parameter is 0
parameter = self._transformed_parameters[orb_index + 1]
if parameter == 0:
continue
# Add the the list of sites and the orbit to the respective cluster maps
cluster_sites = [site.index for site in cluster.lattice_sites]
cluster_to_sites_map.append(cluster_sites)
cluster_to_orbit_map.append(orb_index)
# calculate the number of clusters per orbit
nclusters_per_orbit = [cluster_to_orbit_map.count(i) for i in
range(cluster_to_orbit_map[-1] + 1)]
nclusters_per_orbit = [1] + nclusters_per_orbit
self._cluster_to_sites_map = cluster_to_sites_map
self._cluster_to_orbit_map = cluster_to_orbit_map
self._nclusters_per_orbit = nclusters_per_orbit
def _get_total_energy(self, cluster_instance_activities: List[int]) -> List[float]:
r"""
Calculates the total energy using the expression based on binary variables.
.. math::
H({\boldsymbol x}, {\boldsymbol E})=E_0+
\sum\limits_j\sum\limits_{{\boldsymbol c}
\in{\boldsymbol C}_j}E_jy_{{\boldsymbol c}},
where (:math:`y_{{\boldsymbol c}}=
\prod\limits_{i\in{\boldsymbol c}}x_i`).
Parameters
----------
cluster_instance_activities
list of cluster instance activities, (:math:`y_{{\boldsymbol c}}`)
"""
E = [0.0 for _ in self._transformed_parameters]
for i in range(len(cluster_instance_activities)):
orbit = self._cluster_to_orbit_map[i]
E[orbit + 1] = E[orbit + 1] + cluster_instance_activities[i]
E[0] = 1
E = [0.0 if np.isclose(self._transformed_parameters[orbit], 0.0) else
E[orbit] * self._transformed_parameters[orbit] / self._nclusters_per_orbit[orbit]
for orbit in range(len(self._transformed_parameters))]
return E
[docs]
def get_ground_state(self,
species_count: Dict[str, int] = None,
max_seconds: float = inf,
threads: int = 0) -> Atoms:
"""Finds the ground state for a given structure and species count. If
:attr:`species_count` is not provided when initializing the
instance of this class the first species in the list of
chemical symbols for the active sublattice will be used.
Parameters
----------
species_count
Dictionary with count for one of the species on each active
sublattice. If no count is provided for a sublattice, the
concentration is allowed to vary.
max_seconds
Maximum runtime in seconds.
threads
Number of threads to be used when solving the problem, given that a
positive integer has been provided. If set to :math:`0` the solver default
configuration is used while :math:`-1` corresponds to all available
processing cores.
"""
if species_count is None:
species_count = {}
# Check that the species_count is consistent with the cluster space
all_active_species = set.union(*self._active_species)
for symbol in species_count:
if symbol not in all_active_species:
raise ValueError('The species {} is not present on any of the active sublattices'
' ({})'.format(symbol, self._active_species))
# The model is solved using python-MIPs choice of solver, which is
# Gurobi, if available, and COIN-OR Branch-and-Cut, otherwise.
model = self._model
# Update the species counts
for i, species in enumerate(self._active_species):
count_symbol = self._count_symbols[i]
max_count = len(self._active_sublattices[i].indices)
symbols_to_add = set.intersection(set(species_count), set(species))
if len(symbols_to_add) > 1:
raise ValueError('Provide counts for at most one of the species on each active '
'sublattice ({}), not {}!'.format(self._active_species,
list(species_count)))
elif len(symbols_to_add) == 1:
sym = symbols_to_add.pop()
count = species_count[sym]
if count < 0 or count > max_count:
raise ValueError('The count for species {} ({}) must be a positive integer and'
' cannot exceed the number of sites on the active sublattice '
'({})'.format(sym, count, max_count))
if sym == count_symbol:
xcount = count
else:
xcount = max_count - count
max_slack = 0
else:
xcount = max_slack = max_count
model.constr_by_name('{} count'.format(count_symbol)).rhs = xcount
model.constr_by_name('{} slack'.format(count_symbol)).rhs = max_slack
# Set the number of threads
model.threads = threads
# Optimize the model
self._optimization_status = model.optimize(max_seconds=max_seconds)
# The status of the solution is printed to the screen
if str(self._optimization_status) != 'OptimizationStatus.OPTIMAL':
if str(self._optimization_status) == 'OptimizationStatus.FEASIBLE':
logger.warning('Solution optimality not proven.')
else:
raise Exception('Optimization failed ({0})'.format(str(self._optimization_status)))
# Each of the variables is printed with it's resolved optimum value
gs = self.structure.copy()
active_index_to_sublattice_map = {i: j for j, subl in enumerate(self._active_sublattices)
for i in subl.indices}
for v in model.vars:
if 'atom' in v.name:
index = int(v.name.split('_')[-1])
sublattice_index = active_index_to_sublattice_map[index]
gs[index].symbol = self._reverse_id_maps[sublattice_index][int(v.x)]
# Assert that the solution agrees with the prediction
prediction = self._cluster_expansion.predict(gs)
assert abs(model.objective_value - prediction) < 1e-6
return gs
@property
def optimization_status(self) -> mip.OptimizationStatus:
""" Optimization status. """
return self._optimization_status
@property
def model(self) -> mip.Model:
""" Python-MIP model. """
return self._model.copy()