{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Nanowire Design\n",
"This notebook serves as an example application of the MatOpt framework.\n",
"We consider an example optimization problem of designing a semiconductor compound nanowire."
]
},
{
"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. For semiconductor materials, **WurtziteLattice** is an appropriate type 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 that aligns the resulting lattice with certain orientations expressed by hexagonal Miller Indices. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"IAD = 3.7265 # Angstrom\n",
"orientation = '0001'\n",
"lattice = WurtziteLattice.alignedWith(IAD, orientation)"
]
},
{
"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 linearly periodic, so we will define a **Tiling** object to hold information about the periodicity. In this example, **Cylinder** and **LinearTiling** are the appropriate types, respectively.\n",
"\n",
"Note that we shift the shape of our design space slightly, in order to avoid confusion about which lattice sites (out of those that lie perfectly on the shape facet) should be included."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"nAtomRadius = 6\n",
"nAtomUnitLength = 2\n",
"origin = np.zeros(3, dtype=float)\n",
"axisDirection = np.array([0, 0, 1], dtype=float)\n",
"radius = lattice.getShellSpacing(orientation) * (nAtomRadius - 1)\n",
"height = lattice.getLayerSpacing(orientation) * lattice.getUniqueLayerCount(orientation) * nAtomUnitLength\n",
"shape = Cylinder(origin, radius, height, axisDirection)\n",
"shape.shift(-0.001 * shape.Vh)\n",
"tiling = LinearTiling.fromCylindricalShape(shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given the parameters for a design space, we can construct a **Canvas** object to hold information about lattice locations and nearest neighbors. In this example, the object is first constructed from `lattice` and `shape` and then made periodic with `tiling`. A Canvas can be constructed and manipulated via user-defined algorithms."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"canvas = Canvas.fromLatticeAndShape(lattice, shape)\n",
"canvas.makePeriodic(tiling, lattice.getNeighbors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `canvas` object holds information about the design space and the lattice sites, but it does not specify any material building block information. In this example, we are considering an InAs nanowire. To represent material configurations, we will use a **Design** object. We first construct an empty Design from `canvas` and then place **Atom** (i.e., building block) objects into this Design via methods provided by `lattice`.\n",
"\n",
"To debug our work so far, we can create material structure files to load and plot with standard visualization tools. Here, we create a PDB (protein data bank format, www.rcsb.org) file for the Design object. These files can be plotted with visualization packages such as AtomEye or OVITO."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"design = Design(canvas)\n",
"lattice.setDesign(design, Atom('In'), Atom('As'))\n",
"design.toPDB('canvas.pdb')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building a Model\n",
"\n",
"In this example, we will build a model that maximizes the cohesive energy of the nanowire as an indicator of it's thermal stability. To begin, we start by creating a **MatOptModel** object to hold information about the model. Notice that we use a list of empty atoms as our set of building blocks because we will use a cohesive energy function for semiconductor materials that is indepedent of atom types. Thus, we do not need to specify the types of building blocks in our model."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"m = MatOptModel(canvas, [Atom('')])"
]
},
{
"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 fix the core layers of the nanowire specified by the ratio of the core radius and the canvas radius."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"coreRatio = 0.2\n",
"CoreLayers = [i for i, p in enumerate(canvas.Points) if p[0] ** 2 + p[1] ** 2 < (coreRatio * radius) ** 2]\n",
"m.Yi.rules.append(FixedTo(1, sites=CoreLayers))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we introduce constraints that require atoms to be placed on top of each other radially, avoiding hollow pockets below the surface. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"CanvasMinusCoreLayers = [i for i, p in enumerate(canvas.Points) if i not in CoreLayers]\n",
"NeighborsInside = [[j for j in canvas.NeighborhoodIndexes[i] if (\n",
" j is not None and canvas.Points[j][0] ** 2 + canvas.Points[j][1] ** 2 <\n",
" p[0] ** 2 + p[1] ** 2 - DBL_TOL)] for i, p in enumerate(canvas.Points)]\n",
"m.Yi.rules.append(ImpliesNeighbors(concs=(m.Yi, GreaterThan(1)),\n",
" sites=CanvasMinusCoreLayers,\n",
" neighborhoods=NeighborsInside))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we define a descriptor for each site's contribution to the overall cohesion as a function of each site's coordination environment. We derive this individual atom contribution from the Khor-Das Sarma-type empirical potential energy function. The parameters come from the literature. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"CNBounds = (0, 4)\n",
"p = -0.22479870084561238\n",
"q = -0.9092660150083058\n",
"alpha = -0.3684513\n",
"BPs = list(range(CNBounds[0], CNBounds[1] + 1))\n",
"Vals = [(p * pow(cn, 1 - alpha) - q * cn) for cn in BPs]\n",
"m.addSitesDescriptor('Vi', bounds=(min(Vals), max(Vals)),\n",
" rules=PiecewiseLinear(values=Vals, breakpoints=BPs, input_desc=m.Ci, con_type='UB'))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we define the cohesive energy as a sum of contributions from all `Vi` descriptors."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"size = 216\n",
"m.addGlobalDescriptor('Ecoh', rules=EqualTo(SumSites(desc=m.Vi,coefs=(1.0 / size))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we add constraints on the size and composition of the resulting designs."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"m.addGlobalDescriptor('Size', bounds=(size, size), rules=EqualTo(SumSites(desc=m.Yi)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving the Model\n",
"\n",
"Given a fully formed Pyomo model, we have several capabilities to optimize and visualize the solution. \n",
"In this example, we simply call the `maximize` method to optimize the cohesive energy."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer 12.10.0.0\n",
" with Simplex, Mixed Integer & Barrier Optimizers\n",
"5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21\n",
"Copyright IBM Corp. 1988, 2019. All Rights Reserved.\n",
"\n",
"Type 'help' for a list of available commands.\n",
"Type 'help' followed by a command name for more\n",
"information on commands.\n",
"\n",
"CPLEX> Logfile 'cplex.log' closed.\n",
"Logfile 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmpzlk98pcl.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\\tmp6kc2i4wc.pyomo.lp' read.\n",
"Read time = 0.02 sec. (0.79 ticks)\n",
"CPLEX> Problem name : C:\\Users\\xiang\\AppData\\Local\\Temp\\tmp6kc2i4wc.pyomo.lp\n",
"Objective sense : Maximize\n",
"Variables : 3439 [Nneg: 1, Fix: 1, Box: 520, Free: 1,\n",
" Binary: 2400, General Integer: 516]\n",
"Objective nonzeros : 1\n",
"Linear constraints : 9119 [Less: 6516, Greater: 2080, Equal: 523]\n",
" Nonzeros : 22127\n",
" RHS nonzeros : 3498\n",
"\n",
"Variables : Min LB: 0.000000 Max UB: 216.0000 \n",
"Objective nonzeros : Min : 1.000000 Max : 1.000000 \n",
"Linear constraints :\n",
" Nonzeros : Min : 0.004629630 Max : 1.000000 \n",
" RHS nonzeros : Min : 0.1308177 Max : 16.00000 \n",
"CPLEX> Version identifier: 12.10.0.0 | 2019-11-26 | 843d4de2ae\n",
"CPXPARAM_TimeLimit 360\n",
"CPXPARAM_MIP_Tolerances_AbsMIPGap 0\n",
"CPXPARAM_MIP_Tolerances_MIPGap 0\n",
"Tried aggregator 2 times.\n",
"MIP Presolve eliminated 4306 rows and 1351 columns.\n",
"MIP Presolve modified 1344 coefficients.\n",
"Aggregator did 48 substitutions.\n",
"Reduced MIP has 4765 rows, 2040 columns, and 11724 nonzeros.\n",
"Reduced MIP has 1008 binaries, 516 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.03 sec. (22.50 ticks)\n",
"Found incumbent of value 0.108386 after 0.09 sec. (47.27 ticks)\n",
"Probing time = 0.03 sec. (9.51 ticks)\n",
"Tried aggregator 1 time.\n",
"Detecting symmetries...\n",
"MIP Presolve eliminated 171 rows and 57 columns.\n",
"Reduced MIP has 4594 rows, 1983 columns, and 11325 nonzeros.\n",
"Reduced MIP has 951 binaries, 516 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.02 sec. (12.50 ticks)\n",
"Probing time = 0.03 sec. (8.78 ticks)\n",
"Clique table members: 7696.\n",
"MIP emphasis: balance optimality and feasibility.\n",
"MIP search method: dynamic search.\n",
"Parallel mode: deterministic, using up to 8 threads.\n",
"Root relaxation solution time = 0.23 sec. (194.34 ticks)\n",
"\n",
" Nodes Cuts/\n",
" Node Left Objective IInf Best Integer Best Bound ItCnt Gap\n",
"\n",
"* 0+ 0 0.1084 4.8576 --- \n",
" 0 0 2.3468 1363 0.1084 2.3468 11 --- \n",
"* 0+ 0 1.7329 2.3468 35.43%\n",
"* 0+ 0 1.8527 2.3468 26.67%\n",
" 0 0 2.1194 1363 1.8527 Cuts: 1248 1600 14.40%\n",
"* 0+ 0 1.9802 2.1194 7.03%\n",
" 0 0 2.0221 1363 1.9802 Cuts: 775 2508 2.12%\n",
" 0 0 2.0181 1363 1.9802 Cuts: 34 2800 1.92%\n",
" 0 0 2.0176 1363 1.9802 Cuts: 32 2873 1.89%\n",
" 0 0 2.0174 1363 1.9802 Cuts: 38 2932 1.88%\n",
"Detecting symmetries...\n",
" 0 0 2.0174 1363 1.9802 Cuts: 11 2943 1.88%\n",
"* 0+ 0 1.9823 2.0174 1.77%\n",
"* 0+ 0 1.9979 2.0174 0.97%\n",
"Detecting symmetries...\n",
" 0 2 2.0174 1252 1.9979 2.0174 2943 0.97%\n",
"Elapsed time = 3.16 sec. (2061.34 ticks, tree = 0.02 MB, solutions = 6)\n",
" 8 5 2.0007 897 1.9979 2.0172 3749 0.97%\n",
"\n",
"Clique cuts applied: 4\n",
"Implied bound cuts applied: 562\n",
"Flow cuts applied: 174\n",
"Mixed integer rounding cuts applied: 72\n",
"Zero-half cuts applied: 19\n",
"Lift and project cuts applied: 1\n",
"Gomory fractional cuts applied: 1\n",
"\n",
"Root node processing (before b&c):\n",
" Real time = 3.13 sec. (2048.62 ticks)\n",
"Parallel b&c, 8 threads:\n",
" Real time = 0.69 sec. (458.60 ticks)\n",
" Sync time (average) = 0.35 sec.\n",
" Wait time (average) = 0.00 sec.\n",
" ------------\n",
"Total (root+branch&cut) = 3.81 sec. (2507.22 ticks)\n",
"\n",
"Solution pool: 6 solutions saved.\n",
"\n",
"MIP - Integer optimal solution: Objective = 1.9979448349e+00\n",
"Solution time = 3.81 sec. Iterations = 9300 Nodes = 66\n",
"Deterministic time = 2507.23 ticks (657.55 ticks/sec)\n",
"\n",
"CPLEX> Incumbent solution written to file 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmpl2kim4rh.cplex.sol'.\n",
"CPLEX> The solver exited normally.\n",
"A feasible and provably optimal solution is available.\n",
"The Design has objective: 1.9979448349335396\n"
]
}
],
"source": [
"optimalDesign = None\n",
"try:\n",
" optimalDesign = m.maximize(m.Ecoh, tilim=360, solver='cplex')\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. In this case, we will first label atoms with their identities and then replicate the design four times to better visualize the periodic pattern. "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"if optimalDesign is not None:\n",
" for i, p in enumerate(optimalDesign.Canvas.Points):\n",
" if optimalDesign.Contents[i] is not None:\n",
" if lattice.isASite(p):\n",
" optimalDesign.setContent(i, Atom('In'))\n",
" elif lattice.isBSite(p):\n",
" optimalDesign.setContent(i, Atom('As'))\n",
" optimalDesign.toPDB('result.pdb')\n",
" periodicDesign = deepcopy(optimalDesign)\n",
" for k in range(4):\n",
" for i, p in enumerate(optimalDesign.Canvas.Points):\n",
" periodicDesign.add(p + (k + 1) * shape.Vh, optimalDesign.Contents[i])\n",
" for k in range(4):\n",
" for i, p in enumerate(optimalDesign.Canvas.Points):\n",
" periodicDesign.add(p - (k + 1) * shape.Vh, optimalDesign.Contents[i])\n",
" periodicDesign.toPDB('periodic_result.pdb')"
]
}
],
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}