Bifunctional Surface Design¶
Bifunctional Catalyst Design¶
This notebook serves as an example application of the MatOpt framework. We consider an example optimization problem of designing a nanostructured bifunctional catalyst. This example is a simplified representation of the system presented in [1].
[1] Nunez, M., & Vlachos, D. G. (2019). Ind. Eng. Chem. Res., 58, 6146-6154.
Importing Packages¶
We start by importing several standard Python modules for convienience.
from copy import deepcopy
import numpy as np
Finally, we import the MatOpt package in its entirety.
from idaes.apps.matopt import *
Representing Materials¶
To begin, we define a Lattice object. In this example, FCCLattice is a child class of Lattice. This object will serve to define neighbor connections and helps us generically create other objects.
IAD = 2.828 # Angstrom
Lat = FCCLattice.alignedWith111(IAD)
Next, we define a Shape object that we will use to specify a design space. Additionally, in this example our design space is periodic, so we will define a Tiling object to hold information about the periodicity. In this example, Parallelepiped and PlanarTiling are the appropriate child classes for these objects, respectively.
Note that we shift the shape of our design space slightly, in order to avoid confusion about which lattice sites that lie perfectly on the shape facet should be included.
nUnitCellsOnEdge = 8
nLayers = 4
a = nUnitCellsOnEdge*IAD
b = a
c = nLayers*Lat.FCC111LayerSpacing
alpha = np.pi/2
beta = np.pi/2
gamma = np.pi/3
S = Parallelepiped.fromEdgesAndAngles(a,b,c,alpha,beta,gamma)
S.shift(np.array([-0.01*a,-0.01*b,-0.01*c]))
T = PlanarTiling(S)
Given the parameters for a design space, we can construct a Canvas object to hold information about points and nearest neighbors. In this example, the object is efficiently constructed from a scan over lattice sites. In general, the Canvas can be constructed and manipulated via user-defined algorithms.
Canv = Canvas.fromLatticeAndTilingScan(Lat,T)
The Canvas object hold information about the design space and the lattice sites, but it does not specify any material building block information. To represent material configurations, use a Design object.
Initially, the Design is empty. There are several ways to place Atom (i.e., building block) objects in a Design. In this example, we are initialize the Design to hold all Pt atoms.
To debug our work so far, we can create material structure files to load and plot with standard visualization tools such as AtomEye. Here, we create PDB (protein data bank format, www.rcsb.org) and CFG (AtomEye configuration, li.mit.edu/A/Graphics/A/) files for the undoped design.
D = Design(Canv,Atom('Pt'))
D.toPDB('canvas.pdb')
D.toCFG('canvas.cfg',BBox=S)
Representing Conformations¶
In this material system, we would like to model the presence of facet and edge sites on a patchy bimetallic catalyst surface. To do this generically, we will create a list of conformations. This list will later be used by MatOpt modeling methods to create common descriptor formulations.
To begin, we create another Canvas object with one shell of neighbors around a lattice location. Then, we create a list of Designs and set their contents to match our intended conformations. To debug our work, we also output conformations to file for plotting.
MotifCanvas = Canvas()
MotifCanvas.addLocation(np.array([0,0,0],dtype=float),NNeighbors=12)
MotifCanvas.addShell(Lat.getNeighbors)
Confs = [[None]*len(MotifCanvas.NeighborhoodIndexes[0]) for _ in range(7)]
iToSetNi = [[3,4,5,6,7,8],
[3,4,5,6],
[4,5,6,7],
[5,6,7,8],
[6,7,8,3],
[7,8,3,4],
[8,3,4,5]]
iToSetPt = [[9,10,11],
[9,10,11],
[9,10,11],
[9,10,11],
[9,10,11],
[9,10,11],
[9,10,11]]
for iConf,Conf in enumerate(Confs):
for i in iToSetNi[iConf]:
Conf[i] = Atom('Ni')
for i in iToSetPt[iConf]:
Conf[i] = Atom('Pt')
Building the Model¶
To begin, we define several sets and constants that will be used in creating the model.
TypeAConfs = [0]
TypeBConfs = [1,2,3,4,5,6]
LocsToFixPt = [i for i in range(len(Canv)) if Canv.Points[i][2] < Lat.FCC111LayerSpacing*2.5]
LocsToExcludePt = [i for i in range(len(Canv)) if i not in LocsToFixPt]
CanvTwoBotLayers = [i for i in range(len(Canv)) if Canv.Points[i][2] < Lat.FCC111LayerSpacing*1.5]
CanvMinusTwoBotLayers = [i for i in range(len(Canv)) if i not in CanvTwoBotLayers]
OneLocToFix = [min(LocsToExcludePt)]
TileSizeSquared = nUnitCellsOnEdge**2
CatNorm = TileSizeSquared*6.0
UndefectedSurfE = 0.129758
maxSurfE = 999
CatWeight = 1.0
Atoms = [Atom('Ni'),Atom('Pt')]
Next, we create a *MatOptModel* object.
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.
First, we fix the composition of atoms in the appropriate layers. Effectively, we are designing the defects in a single layer of Ni on top of an undefected Pt surface.
m.Yik.rules.append(FixedTo(1,sites=LocsToFixPt,site_types=[Atom('Pt')]))
m.Yik.rules.append(FixedTo(0,sites=LocsToExcludePt,site_types=[Atom('Pt')]))
Next, we define indicators for the presence of groups of conformations (corresponding to facet and edge sites) in the design. We arbitrarily fix one site to be a facet-type site, breaking symmetry and improving the tractability of the resulting optimization models.
m.Zic.rules.append(FixedTo(1,sites=OneLocToFix,confs=TypeAConfs))
m.Zic.rules.append(Implies(concs=(m.Yik,EqualTo(1,site_types=[Atom('Ni')]))))
SumAConfsExpr = SumConfs(m.Zic,confs_to_sum=TypeAConfs)
SumBConfsExpr = SumConfs(m.Zic,confs_to_sum=TypeBConfs)
m.addBondsDescriptor('SiteCombinations',binary=True,
rules=ImpliesSiteCombination(Canv,
(SumAConfsExpr,GreaterThan(1)),
(SumBConfsExpr,GreaterThan(1))))
Next, we define activity as a normalized sum of contributions from site combinations. Additionally, we introduce a model for the surface energy of sites as a piecewise linear function of coordination number.
m.addGlobalDescriptor('Activity',
rules=EqualTo(SumBonds(m.SiteCombinations,coefs=1/CatNorm)))
EiVals = [0, -0.04293*3+0.41492, -0.04293*10+0.41492, 0.05179*11-0.62148, 0]
EiBPs = [0, 3, 10, 11, 12]
m.addSitesDescriptor('Ei',
rules=PiecewiseLinear(values=EiVals,
breakpoints=EiBPs,
input_desc=m.Ci),
sites=CanvMinusTwoBotLayers)
m.addGlobalDescriptor('Esurf',
rules=EqualTo(SumSites(m.Ei,coefs=1/TileSizeSquared,offset=0.101208)))
m.addGlobalDescriptor('Stability',
rules=EqualTo(LinearExpr(m.Esurf,1/UndefectedSurfE)))
Finally, we introduce a single descriptor for the weighted combination of acitivity and stability. By changing the parameter weighting the catalytic portion of the objective function, we can optimize for a range of designs optimizing stability and activity.
m.addGlobalDescriptor('ActAndStab',
rules=EqualTo(LinearExpr(descs=[m.Stability,m.Activity],
coefs=[-(1-CatWeight),CatWeight])))
Solving the Model¶
Given a fully formed Pyomo model, we have several capabilities to optimize and visualize the solution. In this example, we simply call the maximize method to optimize the balance of activity and stability
D = None
try:
D = m.maximize(m.ActAndStab,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/tmpu_dkga8c.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/tmp62lecdrf.pyomo.lp' read.
Read time = 0.06 sec. (3.48 ticks)
CPLEX> Problem name : /tmp/tmp62lecdrf.pyomo.lp
Objective sense : Maximize
Variables : 8194 [Nneg: 1, Free: 644, Binary: 7165,
General Integer: 384]
Objective nonzeros : 1
Linear constraints : 36613 [Less: 35200, Greater: 128, Equal: 1285]
Nonzeros : 92491
RHS nonzeros : 12680
Variables : Min LB: 0.000000 Max UB: 12.00000
Objective nonzeros : Min : 1.000000 Max : 1.000000
Linear constraints :
Nonzeros : Min : 0.002604167 Max : 12.00000
RHS nonzeros : Min : 0.1012080 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¶
Once the model is solved, we can interpret the solutions as labelings of a Design object. To accompolish this, we use the *setDesignFromModel* function. Then, we can write the Design object to PDB or CFG files for plotting.
if(D is not None):
D.toCFG('result.cfg',BBox=S)
PeriodicD = T.replicateDesign(D,4)
PeriodicS = deepcopy(S)
PeriodicS.scale(np.array([4,4,1]))
PeriodicD.toCFG('periodic_result.cfg',BBox=PeriodicS)