Creating a basic Cluster Expansion#
[1]:
import random
import numpy as np
from monty.serialization import loadfn, dumpfn
from pymatgen.core.structure import Structure
from smol.cofe import ClusterSubspace, StructureWrangler, ClusterExpansion, RegressionData
[2]:
# load the prim structure
lno_prim = loadfn('data/lno_prim.json')
# load the fitting data
lno_entries = loadfn("data/lno_entries.json")
0) The prim structure#
The prim structure defines the configurational space for the Cluster Expansion. The configurational space is defined by the site compositional spaces and the crystal symetries of the prim structure. The occupancy of the sites determine site compositional spaces. Sites are active if they have compositional degrees of freedom.
Active sites have fractional compositions. Vacancies are allowed in sites where the composition does not sum to one.
Is active. The allowed species are: Li+ and vacancies.
Is active. The allowed species are: Ni3+ and Ni4+.
Is not active. Only O2- is allowed.
Is not active. Only O2- is allowed.
[3]:
print(lno_prim)
Full Formula (Li0.5 Ni1 O2)
Reduced Formula: Li0.5Ni1O2
abc : 2.969848 2.969848 5.143928
angles: 73.221350 73.221347 60.000002
pbc : True True True
Sites (4)
# SP a b c
--- ---------------------- ---- ---- ----
0 Li+:0.500 0 0 0
1 Ni3+:0.500, Ni4+:0.500 0.5 0.5 0.5
2 O2- 0.75 0.75 0.75
3 O2- 0.25 0.25 0.25
1) The cluster subspace#
The ClusterSubspace represents all the orbits (groups of equivalent clusters) that will be considered when fitting the cluster expansion. Its main purpose is to compute the correlations functions for each included orbit given a structure in the compositional space defined by the prim.
In order to do be able to compute the correlation functions, the given structure must match the prim structure in a “crystallographic” sense but allowing for compositional degrees of freedom in the “active” sites.
A cluster subspace most easily created by providing: 1. The prim structure representing the configurational space. 2. A set of diameter cutoffs for each size of orbit we want to consider. 3. A type of site basis function to use.
There are more options allowed by the code to fine grain and tune. See other notebooks for advanced use cases.
[4]:
subspace = ClusterSubspace.from_cutoffs(
lno_prim,
cutoffs={2: 5, 3: 4.1}, # will include orbits of 2 and 3 sites.
basis='sinusoid', # sets the site basis type, default is indicator
supercell_size='O2-'
)
# supercell_size specifies the method to determine the supercell size
# when trying to match a structure.
# (See pymatgen.structure_matcher.StructureMatcher for more info)
print(subspace) # single site and empty orbits are always included.
Basis/Orthogonal/Orthonormal : sinusoid/True/True
Unit Cell Composition : Li+0.5 Ni3+0.5 Ni4+0.5 O2-2
Number of Orbits : 11
No. of Correlation Functions : 11
Cluster Cutoffs : 2: 4.20, 3: 2.97
External Terms : []
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 2.9698 6 1 |
| 4 2 2.9698 3 1 |
| 5 2 2.9698 3 1 |
| 6 2 4.2000 6 1 |
| 7 3 2.9698 6 1 |
| 8 3 2.9698 6 1 |
| 9 3 2.9698 2 1 |
| 10 3 2.9698 2 1 |
------------------------------------------------------------------------
1.1) Computing a correlation vector.#
A correlation vector for a specific structure (represents the feature vector) used to train and predict target values.
[5]:
structure = lno_entries[1].structure
corr = subspace.corr_from_structure(structure)
print(f'The correlation vector for a structure with'
f' composition {structure.composition} is: '
f'\n{corr}')
The correlation vector for a structure with composition Li+1 Ni4+5 Ni3+1 O2-12 is:
[1. 0.66666667 0.66666667 0.44444444 0.55555556 0.55555556
0.55555556 0.44444444 0.44444444 0.66666667 0.66666667]
2) The structure wrangler#
The StructureWrangler is a class that will is used to create and organize the data that will be used to train (and possibly test) the cluster expansion. It makes sure that all the supplied structures appropriately match the prim structure, and obtains the necessary information to correctly normalize target properties (such as energy) necessary for training.
Training data is added to a StructureWrangler using ComputedStructureEntry instances from pymatgen.
Matching relaxed structures can be a tricky problem, especially for ionic systems with vacancies. See the notebook on structure matching for tips on how to tweak parameters.
[6]:
wrangler = StructureWrangler(subspace)
# the energy is taken directly from the ComputedStructureEntry
# any additional properties can also be added, see notebook on
# training data preparation for an example.
for entry in lno_entries:
wrangler.add_entry(entry, verbose=True)
# The verbose flag will print structures that fail to match.
print(f'\nTotal structures that match {wrangler.num_structures}/{len(lno_entries)}')
/home/lbluque/Develop/smol/smol/cofe/wrangling/wrangler.py:770: UserWarning: Unable to match Ni4+6 O2-12 with energy -188.28833 to supercell_structure. Throwing out.
Error Message: Supercell could not be found from structure
warnings.warn(
/home/lbluque/Develop/smol/smol/cofe/wrangling/wrangler.py:770: UserWarning: Unable to match Li+2 Ni4+4 Ni3+2 O2-12 with energy -200.13866 to supercell_structure. Throwing out.
Error Message: Mapping could not be found from structure.
warnings.warn(
/home/lbluque/Develop/smol/smol/cofe/wrangling/wrangler.py:770: UserWarning: Unable to match Li+2 Ni3+2 Ni4+4 O2-12 with energy -200.42049 to supercell_structure. Throwing out.
Error Message: Mapping could not be found from structure.
warnings.warn(
/home/lbluque/Develop/smol/smol/cofe/wrangling/wrangler.py:770: UserWarning: Unable to match Li+3 Ni4+4 Ni2+1 Ni3+1 O2-12 with energy -206.70884 to supercell_structure. Throwing out.
Error Message: Supercell could not be found from structure
warnings.warn(
Total structures that match 27/31
Training a cluster expansion is one of the most critical steps. This is how you get effective cluster interactions (ECI’s). To do so you need an estimator class that implements some form of regression model. In this case we will use simple least squares regression using the LinearRegression estimator from scikit-learn.
In smol the coefficients from the fit are not exactly the ECI’s but the ECI times the multiplicity of their orbit.
[7]:
from sklearn.linear_model import LinearRegression
# Set fit_intercept to False because we already do this using
# the empty cluster.
estimator = LinearRegression(fit_intercept=False)
estimator.fit(wrangler.feature_matrix, wrangler.get_property_vector('energy'))
coefs = estimator.coef_
3.1) Check the quality of the fit#
There are many ways to evaluate the quality of a fit. The simplest involve stadard training set prediction error metrics. But when evaluating a CE more seriously we need to consider further metrics and how the CE will be used. Here we will just look at in sample mean squared error and max error.
[8]:
from sklearn.metrics import mean_squared_error, max_error
train_predictions = np.dot(wrangler.feature_matrix, coefs)
rmse = mean_squared_error(
wrangler.get_property_vector('energy'), train_predictions, squared=False
)
maxer = max_error(wrangler.get_property_vector('energy'), train_predictions)
print(f'RMSE {1E3 * rmse} meV/prim')
print(f'MAX {1E3 * maxer} meV/prim')
RMSE 11.03007690876058 meV/prim
MAX 22.826186581717423 meV/prim
4) The cluster expansion#
Now we can use the above work to create the ClusterExpansion. The cluster expansion can be used to predict the fitted property for new structures, either for testing quality or for simulations such as in Monte Carlo. Note that when using the predict function, the cluster expansion will have to match the given structure if it has not seen it before. We will also store the details of the regression model used to fit the cluster expansion by using a RegressionData object.
[9]:
reg_data = RegressionData.from_sklearn(
estimator, wrangler.feature_matrix,
wrangler.get_property_vector('energy')
)
expansion = ClusterExpansion(
subspace, coefficients=coefs, regression_data=reg_data
)
structure = random.choice(wrangler.structures)
prediction = expansion.predict(structure, normalized=True)
print(
f'The predicted energy for a structure with composition '
f'{structure.composition} is {prediction} eV/prim.\n'
)
print(f'The fitted coefficients are:\n{expansion.coefs}\n')
print(f'The effective cluster interactions are:\n{expansion.eci}\n')
print(expansion)
The predicted energy for a structure with composition Li+4 Ni3+4 Ni4+2 O2-12 is -35.48913608767167 eV/prim.
The fitted coefficients are:
[-3.44424307e+01 1.52944807e+00 1.52944807e+00 -7.11937730e-02
1.45252212e-01 4.23347433e-02 -9.28828072e-02 1.51736904e-02
-5.89723850e-02 2.69095444e-02 1.10210719e-02]
The effective cluster interactions are:
[-3.44424307e+01 1.52944807e+00 1.52944807e+00 -1.18656288e-02
4.84174038e-02 1.41115811e-02 -1.54804679e-02 2.52894839e-03
-9.82873083e-03 1.34547722e-02 5.51053597e-03]
Basis/Orthogonal/Orthonormal : sinusoid/True/True
Unit Cell Composition : Li+0.5 Ni3+0.5 Ni4+0.5 O2-2
Number of Orbits : 11
No. of Correlation Functions : 11
Cluster Cutoffs : 2: 4.20, 3: 2.97
External Terms : []
Regression Data : estimator=LinearRegression
module=sklearn.linear_model._base
parameters={'copy_X': True, 'fit_intercept': False, 'n_jobs': None, 'positive': False}
Target Property : mean=-34.6449 std=1.2993
ECI-based Property : mean=-34.4424 std=2.1655
Fit Summary
----------------------------------------------------------------------------------------------------
| ID Orbit ID Degree Cluster Diameter ECI Feature AVG Feature STD ECI * STD |
| 0 0 0 NA -34.442 1.000 0.000 -0.000 |
| 1 1 1 0.0000 1.529 -0.062 0.426 0.651 |
| 2 2 1 0.0000 1.529 -0.062 0.426 0.651 |
| 3 3 2 2.9698 -0.012 0.230 0.310 -0.004 |
| 4 4 2 2.9698 0.048 0.218 0.442 0.021 |
| 5 5 2 2.9698 0.014 0.111 0.342 0.005 |
| 6 6 2 4.2000 -0.015 0.337 0.325 -0.005 |
| 7 7 3 2.9698 0.003 -0.070 0.357 0.001 |
| 8 8 3 2.9698 -0.010 -0.029 0.350 -0.003 |
| 9 9 3 2.9698 0.013 -0.062 0.385 0.005 |
| 10 10 3 2.9698 0.006 -0.111 0.416 0.002 |
----------------------------------------------------------------------------------------------------
/home/lbluque/opt/miniconda3/envs/matx_dev/lib/python3.10/site-packages/pymatgen/core/periodic_table.py:1051: UserWarning: Use of properties is now deprecated. Set the spin by setting the spin arg instead.
warnings.warn("Use of properties is now deprecated. Set the spin by setting the spin arg instead.")
5) Saving your work#
All core classes in smol are MSONables and so can be saved using their as_dict methods or better yet with monty.serialization.dumpfn.
Currently there is also a convenience function in smol that will nicely save all of your work for you in a standardized way. Work saved with the save_work function is saved as a dictionary with standardized names for the classes. Since a work flow should only contain 1 of each core classes the function will complain if you give it two of the same class (i.e. two wranglers)
[10]:
from smol.io import save_work
file_path = 'data/basic_ce.mson'
# we can save the subspace as well, but since both the wrangler
# and the expansion have it, there is no need to do so.
save_work(file_path, wrangler, expansion)
5.1) Loading previously saved work#
[11]:
from smol.io import load_work
work = load_work(file_path)
for name, obj in work.items():
print(f'{name}: {type(obj)}\n')
StructureWrangler: <class 'smol.cofe.wrangling.wrangler.StructureWrangler'>
ClusterExpansion: <class 'smol.cofe.expansion.ClusterExpansion'>