{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Surface Design\n",
"This notebook serves as an example application of the MatOpt framework.\n",
"We consider an example optimization problem of designing a monometallic nanostructured catalyst surface.\n",
"\n",
"For more information, see: Hanselman, Christopher L. and Gounaris, Chrysanthos E. \"A mathematical optimization framework for the design of nanopatterned surfaces.\" *AIChE Journal* 62 (2016): 3250-3263. DOI: [10.1002/aic.15359](https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.15359)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importing Packages\n",
"We start by importing several standard Python modules for convienience. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from copy import deepcopy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we import MatOpt."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from idaes.apps.matopt import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Representing Materials\n",
"To begin, we define a **Lattice** object. In this example, **FCCLattice** is the appropriate a child class of Lattice. This object will serve to define neighbor connections and helps us generically create other objects. We construct our lattice from a class method constructor for FCC lattices aligned with the {111} plane. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"IAD = 2.828427 # Angstrom\n",
"Lat = FCCLattice.alignedWith111(IAD)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"nUnitCellsOnEdge = 4\n",
"nLayers = 6\n",
"a = nUnitCellsOnEdge*IAD\n",
"b = a\n",
"c = nLayers*Lat.FCC111LayerSpacing\n",
"alpha = np.pi/2\n",
"beta = np.pi/2\n",
"gamma = np.pi/3\n",
"S = Parallelepiped.fromEdgesAndAngles(a,b,c,alpha,beta,gamma)\n",
"S.shift(np.array([-0.01*a,-0.01*b,-0.01*c]))\n",
"T = PlanarTiling(S)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"Canv = Canvas.fromLatticeAndTilingScan(Lat,T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. \n",
"\n",
"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 considering a Pt surface for the oxygen reduction reaction. We can initialize a design representing the FCC {111} surface by using a standard constructor for the Design object.\n",
"\n",
"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 design. These files can be plotted with visualization packages such as AtomEye or OVITO."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"D = Design(Canv,Atom('Pt'))\n",
"D.toPDB('undefected.pdb')\n",
"D.toCFG('undefected.cfg',GS=1.0,BBox=S)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building a Model\n",
"\n",
"In this example, we will build a model that maximizes the number of sites that are reactive for the oxygen reduction reaction (ORR). More generally, our model will indicate sites that are within a certain tolerance of a target generalized coordination number (GCN). These target sites can also be constrained to lie within minimum and maximum coordination number to be considered surface sites. \n",
"\n",
"Additionally, we model the surface energy of nanostructured designs. This surface energy can be constrained to be below a threshold and can be included in the objective function. We can parametrically optimize the multi-objective optimization problem by defining a weighting, *CatWeight*, that controls how much weight is given to the catalytic activity term in the objective function. A weighting of 1 corresponds to the optimally active material and a weighting of 0 corresponds to the lowest surface energy design. \n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"Atoms = [Atom('Pt')]\n",
"TargetGCN = 8.0\n",
"CNsurfMin = 3\n",
"CNsurfMax = 9\n",
"TileSizeSquared = nUnitCellsOnEdge**2\n",
"UndefectedSurfE = 0.129758\n",
"maxSurfE = 999\n",
"CatWeight = 1.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To begin, we start by creating a ***MatOptModel*** object to hold information about the model. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"m = MatOptModel(Canv,Atoms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we introduce two rules to fix special sites in the design. \n",
"We fix the bottom two layers of atoms to exist, creating underlying bulk layers above which we will introduce nanostruced defets.\n",
"We also fix an arbitrary atom in the top layer, breaking symetry of the design space and resulting in easier to solve opitmization problems without actually restricting the designs that can be possibly represented. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"CanvTwoBotLayers = [i for i in range(len(Canv)) \n",
" if Canv.Points[i][2] < 1.5*Lat.FCC111LayerSpacing]\n",
"CanvMinusTwoBotLayers = [i for i in range(len(Canv)) \n",
" if i not in CanvTwoBotLayers]\n",
"OneSiteInTopLayer = [min([i for i in range(len(Canv)) \n",
" if Canv.Points[i][2] > (nLayers-1.5)*Lat.FCC111LayerSpacing])]\n",
"m.Yi.rules.append(FixedTo(1,sites=OneSiteInTopLayer))\n",
"m.Yi.rules.append(FixedTo(1,sites=CanvTwoBotLayers))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we introduce constraints thtat require atoms to be placed on top of each other, avoiding hollow pockets below the surface. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"NeighborsBelow = [[j for j in Canv.NeighborhoodIndexes[i] \n",
" if(j is not None and\n",
" Canv.Points[j][2] Logfile 'cplex.log' closed.\n",
"Logfile 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmpl2xbudfz.cplex.log' open.\n",
"CPLEX> New value for absolute mixed integer optimality gap tolerance: 0\n",
"CPLEX> New value for mixed integer optimality gap tolerance: 0\n",
"CPLEX> New value for time limit in seconds: 360\n",
"CPLEX> Problem 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmp8in2pgv2.pyomo.lp' read.\n",
"Read time = 0.00 sec. (0.03 ticks)\n",
"CPLEX> Problem name : C:\\Users\\xiang\\AppData\\Local\\Temp\\tmp8in2pgv2.pyomo.lp\n",
"Objective sense : Maximize\n",
"Variables : 184 [Nneg: 1, Box: 9, Free: 42, Binary: 119,\n",
" General Integer: 12, Other: 1]\n",
"Objective nonzeros : 1\n",
"Linear constraints : 469 [Less: 412, Greater: 8, Equal: 49]\n",
" Nonzeros : 960\n",
" RHS nonzeros : 117\n",
"\n",
"Variables : Min LB: 0.0000000 Max UB: 999.0000 \n",
"Objective nonzeros : Min : 1.000000 Max : 1.000000 \n",
"Linear constraints :\n",
" Nonzeros : Min : 0.03741000 Max : 12.00000 \n",
" RHS nonzeros : Min : 0.1012080 Max : 12.00000 \n",
"CPLEX> CPXPARAM_TimeLimit 360\n",
"CPXPARAM_MIP_Tolerances_AbsMIPGap 0\n",
"CPXPARAM_MIP_Tolerances_MIPGap 0\n",
"Tried aggregator 3 times.\n",
"MIP Presolve eliminated 336 rows and 78 columns.\n",
"MIP Presolve modified 51 coefficients.\n",
"Aggregator did 54 substitutions.\n",
"Reduced MIP has 79 rows, 52 columns, and 267 nonzeros.\n",
"Reduced MIP has 29 binaries, 8 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.00 sec. (0.86 ticks)\n",
"Found incumbent of value 0.000000 after 0.00 sec. (1.00 ticks)\n",
"Probing fixed 3 vars, tightened 3 bounds.\n",
"Probing time = 0.00 sec. (0.21 ticks)\n",
"Cover probing fixed 3 vars, tightened 0 bounds.\n",
"Tried aggregator 2 times.\n",
"MIP Presolve eliminated 24 rows and 15 columns.\n",
"MIP Presolve modified 6 coefficients.\n",
"Aggregator did 1 substitutions.\n",
"Reduced MIP has 54 rows, 36 columns, and 161 nonzeros.\n",
"Reduced MIP has 20 binaries, 7 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.00 sec. (0.23 ticks)\n",
"Probing changed sense of 3 constraints.\n",
"Probing time = 0.00 sec. (0.08 ticks)\n",
"Tried aggregator 2 times.\n",
"Aggregator did 3 substitutions.\n",
"Reduced MIP has 51 rows, 33 columns, and 158 nonzeros.\n",
"Reduced MIP has 17 binaries, 7 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.00 sec. (0.18 ticks)\n",
"Probing time = 0.00 sec. (0.07 ticks)\n",
"\n",
"Root node processing (before b&c):\n",
" Real time = 0.00 sec. (1.92 ticks)\n",
"Parallel b&c, 8 threads:\n",
" Real time = 0.00 sec. (0.00 ticks)\n",
" Sync time (average) = 0.00 sec.\n",
" Wait time (average) = 0.00 sec.\n",
" ------------\n",
"Total (root+branch&cut) = 0.00 sec. (1.92 ticks)\n",
"\n",
"Solution pool: 1 solution saved.\n",
"\n",
"MIP - Integer optimal solution: Objective = 2.5000000000e-01\n",
"Solution time = 0.00 sec. Iterations = 0 Nodes = 0\n",
"Deterministic time = 1.92 ticks (1918.70 ticks/sec)\n",
"\n",
"CPLEX> Incumbent solution written to file 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmp09uzd_2e.cplex.sol'.\n",
"CPLEX> The solver exited normally.\n",
"A feasible and provably optimal solution is available.\n",
"The Design has objective: 0.25\n"
]
}
],
"source": [
"D = None\n",
"try:\n",
" D = m.maximize(m.ActAndStab,tilim=360)\n",
"except:\n",
" print('MaOpt can not find usable solver (CPLEX or NEOS-CPLEX)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Processing Solutions\n",
"Once the model is solved, we can plot the resulting design. However, it is often useful to label atoms according to some auxilliary information. In this case, we would like to label atoms that consitute ideal reactive sites. We loop over the sites and set the atom to S to highlight the sites that are reactive. Then, we can write the Design object to PDB or CFG files for plotting.\n",
"\n",
"Additionally, we can manipulate the resulting design to better see the periodic pattern. Here, we replicate the design four times to see the periodic pattern. "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"if(D is not None):\n",
" for i in m.IdealSitei.keys():\n",
" if m.IdealSitei.values[i] > 0.5:\n",
" D.setContent(i,Atom('S'))\n",
" D.toPDB('result.pdb')\n",
" PeriodicD = T.replicateDesign(D,4)\n",
" PeriodicS = deepcopy(S)\n",
" PeriodicS.scale(np.array([4,4,1]))\n",
" PeriodicD.toCFG('periodic_result.cfg',BBox=PeriodicS)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}