Metal Oxide Bulk Design¶
Perovskite Design¶
This notebook serves as an example application of the MatOpt framework on bulk nanostructured materials. We consider the problem of how to optimally place dopant in a perovskite lattice.
For more information, see: Hanselman, Christopher L., et al. "A framework for optimizing oxygen vacancy formation in doped perovskites." Computers & Chemical Engineering 126 (2019): 168-177. DOI: 10.1016/j.compchemeng.2019.03.033
Importing Packages¶
We start by importing MatOpt.
import numpy as np
from idaes.apps.matopt import *
Representing Materials¶
First, we construct a *Lattice* object to hold information about the sites in our material. For this application, we use the *PerovskiteLattice* class with lattice constants A, B, and C.
A = 4.0
B = 4.0
C = 4.0
Lat = PerovskiteLattice(A,B,C)
Next, we construct *Shape* and *Tiling* objects to help define the material locations of interest. In this case, we use *RectPrism* and *CubicTiling*, respectively.
Note that we shift the shape slightly to avoid ambiguity about which sites on the border of the cell should be included.
nUnitCellsOnEdge = 2
S = RectPrism(nUnitCellsOnEdge*A,
nUnitCellsOnEdge*B,
nUnitCellsOnEdge*C)
S.shift(np.array([-0.01,-0.01,-0.01]))
T = CubicTiling(S)
Next, we construct the *Canvas* object which will hold all the information about the sites and information about neighbors. We also define a list of *Atom* objects that serve as the building blocks of our material.
Canv = Canvas.fromLatticeAndTilingScan(Lat,T)
Atoms = [Atom('Ba'),Atom('Fe'),Atom('In'),Atom('O')]
Finally, we load a list of conformations from file that represent a set of dopant configurations that we would like to indicate in the design.
from tempfile import TemporaryDirectory
from zipfile import ZipFile
iDesiredConfs = [394,395,396,397,398,399,400,401,68,69,
70,71,162,163,164,165,166,167,168,169]
with TemporaryDirectory() as ConfDir:
ZipFile('./confs.zip').extractall(ConfDir)
ConfDesigns = loadFromPDBs([str(i) + '.pdb' for i in iDesiredConfs], folder=ConfDir + '/confs/')
Confs = [Conf.Contents for Conf in ConfDesigns]
Building the Model¶
To begin specifying the model, we first define several pieces of information that will help specify the design problem.
Sites = [i for i in range(len(Canv))]
ASites = [i for i in Sites if Lat.isASite(Canv.Points[i])]
BSites = [i for i in Sites if Lat.isBSite(Canv.Points[i])]
OSites = [i for i in Sites if Lat.isOSite(Canv.Points[i])]
pctLocalLB,pctLocalUB = 0,1
pctGlobalLB,pctGlobalUB = 0.0,0.3
LocalBounds = {(i,Atom('In')):(round(pctLocalLB*len(Canv.NeighborhoodIndexes[i])),
round(pctLocalUB*len(Canv.NeighborhoodIndexes[i]))) for i in OSites}
GlobalLB = round(pctGlobalLB*len(BSites))
GlobalUB = round(pctGlobalUB*len(BSites))
Next, we initialize a *MatOptModel* object that will hold all the information about material descriptors and desired functionalities.
m = MatOptModel(Canv,Atoms,Confs)
By default, several basic variables are pre-defined. See the first example, *Monometallic_Nanocluster_Design.ipynb* for a description of basic variables, expressions, and constraint rules.
For this system, we introduce several rules about the allowed placement of atoms in the design. First, we require that all A-sites in the material are occupied by Ba. Next, we require that all O-sites are occupied by O. Thirdly, we forbid Ba and O from being placed in B-sites. And finally, we require that some atom be placed in each B-site. These four rules effectively limit the scope of the optimization to focus on the labeling of B-sites as either Fe or In.
m.Yik.rules.append(FixedTo(1,sites=ASites,site_types=[Atom('Ba')]))
m.Yik.rules.append(FixedTo(1,sites=OSites,site_types=[Atom('O')]))
m.Yik.rules.append(FixedTo(0,sites=BSites,site_types=[Atom('Ba'),Atom('O')]))
m.Yi.rules.append(FixedTo(1,sites=BSites))
To specify additional constraints to the model, we create several descriptors for the activity, local dopant concentration, and the global dopant concentration.
Notice that in each case, we specify a subset of locations or atoms of interest. This is because, for example, our material activity depends on oxygen sites only and it would be nonsensical to try to interpret one of the conformations on a different type of site. Similarly, the dopant budgets are written only over In atoms and not on Ba, Fe, or O.
m.addGlobalDescriptor('Activity',
rules=EqualTo(SumSitesAndConfs(m.Zic,coefs=1/len(OSites),sites_to_sum=OSites)))
m.addGlobalTypesDescriptor('GlobalIndiumConc',bounds=(GlobalLB,GlobalUB),
rules=EqualTo(SumSites(m.Yik,
site_types=[Atom('In')],
sites_to_sum=BSites)))
m.addSitesTypesDescriptor('LocalIndiumConc',bounds=LocalBounds,
rules=EqualTo(SumNeighborSites(m.Yik,
sites=OSites,
site_types=[Atom('In')])))
Solving the Model¶
Given a fully formed model, we can optimize by maximizing or minimizing one of the global descriptors.
D = None
try:
D = m.maximize(m.Activity,tilim=360)
except:
print('MaOpt can not find usable solver (CPLEX or NEOS-CPLEX)')
WARNING: DEPRECATED: SetProduct.set_tuple is deprecated. Use
SetProduct.subsets() to get the operator arguments. (deprecated in 5.7)
(called from /home/runner/.conda/envs/idaes-env/lib/python3.8/site-
packages/idaes/apps/matopt/../matopt/opt/pyomo_modeling.py:324)
Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer Community Edition 12.9.0.0
with Simplex, Mixed Integer & Barrier Optimizers
5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21
Copyright IBM Corp. 1988, 2019. All Rights Reserved.
Type 'help' for a list of available commands.
Type 'help' followed by a command name for more
information on commands.
CPLEX> Logfile 'cplex.log' closed.
Logfile '/tmp/tmp3pd9p547.cplex.log' open.
CPLEX> New value for absolute mixed integer optimality gap tolerance: 0
CPLEX> New value for mixed integer optimality gap tolerance: 0
CPLEX> New value for time limit in seconds: 360
CPLEX> Problem '/tmp/tmp6we5n9qw.pyomo.lp' read.
Read time = 0.01 sec. (0.31 ticks)
CPLEX> Problem name : /tmp/tmp6we5n9qw.pyomo.lp
Objective sense : Maximize
Variables : 843 [Nneg: 1, Box: 25, Free: 1, Binary: 816]
Objective nonzeros : 1
Linear constraints : 1459 [Less: 1424, Equal: 35]
Nonzeros : 17523
RHS nonzeros : 1113
Variables : Min LB: 0.000000 Max UB: 10.00000
Objective nonzeros : Min : 1.000000 Max : 1.000000
Linear constraints :
Nonzeros : Min : 0.04166667 Max : 2.000000
RHS nonzeros : Min : 1.000000 Max : 9.000000
CPLEX> CPLEX Error 1016: Community Edition. Problem size limits exceeded. Purchase at https://ibm.co/2s0wqSa.
Error termination, CPLEX Error 1016.
Solution time = 0.00 sec.
Deterministic time = 0.00 ticks (0.00 ticks/sec)
CPLEX> CPLEX Error 1217: No solution exists.
No file written.
CPLEX> ERROR: evaluating object as numeric value: obj
(object: <class 'pyomo.core.base.objective.ScalarObjective'>)
No value for uninitialized NumericValue object obj
MaOpt can not find usable solver (CPLEX or NEOS-CPLEX)
Processing Solutions¶
If the optimizer was successful in finding an optimal (or just feasible) solution, we can plot the resulting design to any of several standard file formats. However, it is often useful to modify the design to highlight key features. Here, we label all O-sites that constitute one of the desired conformations by replacing the atom with an S.
if(D is not None):
for i,c in m.Zic.keys():
if(m.Zic.values[i,c] > 0.5):
D.setContent(i,Atom('S'))
D.toCFG('result.cfg',BBox=S)