Running periodic ground-state structure solver#
In this example, we illustrate how to create a periodic ground-state solver for a cluster expansion, and how to obtain ground-states for a specified supercell.
This procedure solves for the ground state of a fixed periodic super-cell specified by the supercell matrix used to create the Ensemble class, this means that the solution obtained corresponds only to an upper bound of the global ground state of the corresponding infinitely sized bulk system.
For more details on the nature of the global problem, finding upper and lower bounds please see: https://doi.org/10.1103/PhysRevB.94.134424
cvxpy and a mixed integer programming solver (such as SCIP, GUROBI, CPLEX) are required. Proprietary solvers such as GUROBI and CPLEX are noteably faster than open source alternatives like SCIP. Also note that in very rare circumstances, using SCIP could be risky as SCIP might fail to satisfy a portion of the integer constraints. Further details of their installation can be found in: https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver
Note that the time to solve grows rapidly with the number of species, the number of terms in the cluster expansion, and the super-cell size. If you searching for the ground-state of a large complex system expect very long runtimes.
[1]:
import numpy as np
import crystal_toolkit
from pymatgen.core import Structure, Lattice
from smol.cofe import ClusterSubspace, ClusterExpansion
from smol.cofe.extern import EwaldTerm
from smol.moca import Ensemble
from smol.capp.generate import PeriodicGroundStateSolver
0) Create a Cluster Subspace based on the disordered structure with an Ewald term#
[2]:
# Hypothetical quarternary system
clof = Structure(
Lattice.cubic(2.0),
[{"Ca2+": 0.5, "Li+":0.5}, {"O2-": 0.5, "F-":0.5}],
[[0,0,0], [0.5, 0.5, 0.5]]
)
[3]:
subspace = ClusterSubspace.from_cutoffs(structure=clof, cutoffs={2: 4.0, 3: 3.0})
subspace.add_external_term(EwaldTerm(eta=None)) # Add the external Ewald Term
[4]:
print(subspace)
Basis/Orthogonal/Orthonormal : indicator/False/False
Unit Cell Composition : Ca2+0.5 Li+0.5 O2-0.5 F-0.5
Number of Orbits : 21
No. of Correlation Functions : 21
Cluster Cutoffs : 2: 4.00, 3: 2.83
External Terms : [EwaldTerm(total)]
Orbit Summary
------------------------------------------------------------------------
| ID Degree Cluster Diameter Multiplicity No. Functions |
| 0 0 NA 0 1 |
| 1 1 0.0000 1 1 |
| 2 1 0.0000 1 1 |
| 3 2 1.7321 8 1 |
| 4 2 2.0000 3 1 |
| 5 2 2.0000 3 1 |
| 6 2 2.8284 6 1 |
| 7 2 2.8284 6 1 |
| 8 2 3.3166 24 1 |
| 9 2 3.4641 4 1 |
| 10 2 3.4641 4 1 |
| 11 2 4.0000 3 1 |
| 12 2 4.0000 3 1 |
| 13 3 2.0000 12 1 |
| 14 3 2.0000 12 1 |
| 15 3 2.8284 12 1 |
| 16 3 2.8284 12 1 |
| 17 3 2.8284 12 1 |
| 18 3 2.8284 12 1 |
| 19 3 2.8284 8 1 |
| 20 3 2.8284 8 1 |
------------------------------------------------------------------------
1) Create random ECIs and ClusterExpansion#
[5]:
coefs = np.random.random(size=subspace.num_corr_functions+1)
coefs[0] = -10
coefs[-1] = 0.1
[ ]:
expansion = ClusterExpansion(subspace, coefs)
2) Create a semi-grand canonical ensemble#
By default, the Ensemble object uses ClusterDecompositionProcessor, which greatly reduces the amount of many-body terms used in optimization pseudo-boolean function. Switching to a ClusterExpansionProcessor is allowed but not recommended.
The upper-bound supercell size is defined within the Ensemble. Here, we use a supercell containing 2 primitive units. A semigrand-canonical ensemble is defined as an Ensemble object with chemical potentials specified.
In this example, we used a relatively simple chemical space and a small supercell for a quick demonstration. The Boolean problem size grows exponentially with increasing composition complexity, therefore, the solver might not be able to solve the upper-bound in a large super-cell with high configurational entropy.
[7]:
chempots = {"Li+": 0.0, "Ca2+": 0.0, "O2-":0.0, "F-": 0.0}
sgc_ensemble = Ensemble.from_cluster_expansion(expansion, np.diag([2, 2, 2]), chemical_potentials=chempots)
3) Create a solver instance from the ensemble#
Charge-balance constraints are included by default. If any other constraint is needed, refer to the documentation of GroundStateSolver.
The default solver is “SCIP”. For other solvers supported by cvxpy, see: https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options
[8]:
# Use 1e-6 as a cutoff to cluster terms. Any term with coefficient lower than 1e-6 will not be included into the optimization.
sgc_solver = PeriodicGroundStateSolver(sgc_ensemble, term_coefficients_cutoff=1e-6)
print("Number of variables:", sgc_solver._canonicals.variables.size)
print("Number of auxiliary variables:", sgc_solver._canonicals.auxiliary_variables.size)
print("Number of constraints:", len(sgc_solver._canonicals.constraints))
Number of variables: 32
Number of auxiliary variables: 960
Number of constraints: 3345
4) Solve the problem in semi-grand canonical ensemble#
[9]:
sgc_solver.solve()
[10]:
print("Ground-state energy, un-normalized(eV):", sgc_solver.ground_state_energy)
print("Ground-state structure:", sgc_solver.ground_state_structure)
sgc_solver.ground_state_structure # needs crystal_toolkit for visualization
Ground-state energy, un-normalized(eV): -97.89342026301097
Ground-state structure: Full Formula (Li2 Ca6 O6 F2)
Reduced Formula: LiCa3O3F
abc : 4.000000 4.000000 4.000000
angles: 90.000000 90.000000 90.000000
pbc : True True True
Sites (16)
# SP a b c
--- ---- ---- ---- ----
0 Ca2+ 0 0 0
1 Li+ 0 0 0.5
2 Ca2+ 0 0.5 0
3 Ca2+ 0 0.5 0.5
4 Li+ 0.5 0 0
5 Ca2+ 0.5 0 0.5
6 Ca2+ 0.5 0.5 0
7 Ca2+ 0.5 0.5 0.5
8 F- 0.25 0.25 0.25
9 O2- 0.25 0.25 0.75
10 O2- 0.25 0.75 0.25
11 O2- 0.25 0.75 0.75
12 O2- 0.75 0.25 0.25
13 O2- 0.75 0.25 0.75
14 O2- 0.75 0.75 0.25
15 F- 0.75 0.75 0.75
If you see this text, the Crystal Toolkit Jupyter Lab
extension is not installed. You can install it by running
"pip install crystaltoolkit-extension"
from the same environment you run "jupyter lab".
This only works in Jupyter Lab 3.x or above.
Structure Summary
Lattice
abc : 4.0 4.0 4.0
angles : 90.0 90.0 90.0
volume : 64.0
A : 4.0 0.0 0.0
B : 0.0 4.0 0.0
C : 0.0 0.0 4.0
pbc : True True True
PeriodicSite: Ca2+ (0.0, 0.0, 0.0) [0.0, 0.0, 0.0]
PeriodicSite: Li+ (0.0, 0.0, 2.0) [0.0, 0.0, 0.5]
PeriodicSite: Ca2+ (0.0, 2.0, 0.0) [0.0, 0.5, 0.0]
PeriodicSite: Ca2+ (0.0, 2.0, 2.0) [0.0, 0.5, 0.5]
PeriodicSite: Li+ (2.0, 0.0, 0.0) [0.5, 0.0, 0.0]
PeriodicSite: Ca2+ (2.0, 0.0, 2.0) [0.5, 0.0, 0.5]
PeriodicSite: Ca2+ (2.0, 2.0, 0.0) [0.5, 0.5, 0.0]
PeriodicSite: Ca2+ (2.0, 2.0, 2.0) [0.5, 0.5, 0.5]
PeriodicSite: F- (1.0, 1.0, 1.0) [0.25, 0.25, 0.25]
PeriodicSite: O2- (1.0, 1.0, 3.0) [0.25, 0.25, 0.75]
PeriodicSite: O2- (1.0, 3.0, 1.0) [0.25, 0.75, 0.25]
PeriodicSite: O2- (1.0, 3.0, 3.0) [0.25, 0.75, 0.75]
PeriodicSite: O2- (3.0, 1.0, 1.0) [0.75, 0.25, 0.25]
PeriodicSite: O2- (3.0, 1.0, 3.0) [0.75, 0.25, 0.75]
PeriodicSite: O2- (3.0, 3.0, 1.0) [0.75, 0.75, 0.25]
PeriodicSite: F- (3.0, 3.0, 3.0) [0.75, 0.75, 0.75]
5) Create and solve a canonical ensemble problem#
Our solver also supports solving the ground-state in canonical ensembles. In doing so, one only needs to create a canonical ensemble (i.e., an ensemble with no chemical potentials provided). When using a canonical ensemble, note that either a fixed composition or an intial occupancy to determine the composition to fix must be provided as an argument to initialize the solver.
[11]:
# Creating the ensemble.
canonical_ensemble = Ensemble.from_cluster_expansion(expansion, np.diag([2, 2, 2]))
for sublattice in canonical_ensemble.sublattices:
print(sublattice)
Sublattice(site_space=Ca2+0.5 Li+0.5 , sites=array([0, 1, 2, 3, 4, 5, 6, 7]), active_sites=array([0, 1, 2, 3, 4, 5, 6, 7]), encoding=array([0, 1]))
Sublattice(site_space=O2-0.5 F-0.5 , sites=array([ 8, 9, 10, 11, 12, 13, 14, 15]), active_sites=array([ 8, 9, 10, 11, 12, 13, 14, 15]), encoding=array([0, 1]))
[12]:
# Fix to Ca4 Li4 O4 F4.
canonical_solver = PeriodicGroundStateSolver(
canonical_ensemble, term_coefficients_cutoff=1e-6,
initial_occupancy=[0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1]
)
print("Number of variables:", canonical_solver._canonicals.variables.size)
print("Number of auxiliary variables:", canonical_solver._canonicals.auxiliary_variables.size)
print("Number of constraints:", len(canonical_solver._canonicals.constraints))
Number of variables: 32
Number of auxiliary variables: 960
Number of constraints: 3348
[13]:
canonical_solver.solve()
[14]:
print("Ground-state energy, un-normalized (eV):", canonical_solver.ground_state_energy)
print("Ground-state structure:", canonical_solver.ground_state_structure)
canonical_solver.ground_state_structure # needs crystal_toolkit for visualization
Ground-state energy, un-normalized (eV): -94.78254966983789
Ground-state structure: Full Formula (Li4 Ca4 O4 F4)
Reduced Formula: LiCaOF
abc : 4.000000 4.000000 4.000000
angles: 90.000000 90.000000 90.000000
pbc : True True True
Sites (16)
# SP a b c
--- ---- ---- ---- ----
0 Li+ 0 0 0
1 Ca2+ 0 0 0.5
2 Ca2+ 0 0.5 0
3 Li+ 0 0.5 0.5
4 Ca2+ 0.5 0 0
5 Li+ 0.5 0 0.5
6 Li+ 0.5 0.5 0
7 Ca2+ 0.5 0.5 0.5
8 O2- 0.25 0.25 0.25
9 F- 0.25 0.25 0.75
10 F- 0.25 0.75 0.25
11 O2- 0.25 0.75 0.75
12 O2- 0.75 0.25 0.25
13 F- 0.75 0.25 0.75
14 F- 0.75 0.75 0.25
15 O2- 0.75 0.75 0.75
If you see this text, the Crystal Toolkit Jupyter Lab
extension is not installed. You can install it by running
"pip install crystaltoolkit-extension"
from the same environment you run "jupyter lab".
This only works in Jupyter Lab 3.x or above.
Structure Summary
Lattice
abc : 4.0 4.0 4.0
angles : 90.0 90.0 90.0
volume : 64.0
A : 4.0 0.0 0.0
B : 0.0 4.0 0.0
C : 0.0 0.0 4.0
pbc : True True True
PeriodicSite: Li+ (0.0, 0.0, 0.0) [0.0, 0.0, 0.0]
PeriodicSite: Ca2+ (0.0, 0.0, 2.0) [0.0, 0.0, 0.5]
PeriodicSite: Ca2+ (0.0, 2.0, 0.0) [0.0, 0.5, 0.0]
PeriodicSite: Li+ (0.0, 2.0, 2.0) [0.0, 0.5, 0.5]
PeriodicSite: Ca2+ (2.0, 0.0, 0.0) [0.5, 0.0, 0.0]
PeriodicSite: Li+ (2.0, 0.0, 2.0) [0.5, 0.0, 0.5]
PeriodicSite: Li+ (2.0, 2.0, 0.0) [0.5, 0.5, 0.0]
PeriodicSite: Ca2+ (2.0, 2.0, 2.0) [0.5, 0.5, 0.5]
PeriodicSite: O2- (1.0, 1.0, 1.0) [0.25, 0.25, 0.25]
PeriodicSite: F- (1.0, 1.0, 3.0) [0.25, 0.25, 0.75]
PeriodicSite: F- (1.0, 3.0, 1.0) [0.25, 0.75, 0.25]
PeriodicSite: O2- (1.0, 3.0, 3.0) [0.25, 0.75, 0.75]
PeriodicSite: O2- (3.0, 1.0, 1.0) [0.75, 0.25, 0.25]
PeriodicSite: F- (3.0, 1.0, 3.0) [0.75, 0.25, 0.75]
PeriodicSite: F- (3.0, 3.0, 1.0) [0.75, 0.75, 0.25]
PeriodicSite: O2- (3.0, 3.0, 3.0) [0.75, 0.75, 0.75]