{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Monometallic Nanocluster Design\n",
"\n",
"In this module, we introduce the **MatOpt** interface for representing material properties and specifying optimization problems. \n",
"\n",
"We have designed the interface with severl goals in mind:\n",
"\n",
"1. To **simplify the representation of nanostructured materials,** streamlining the creation of materials optimization problems. \n",
"2. To provide a simple interface so that users **do not need to understand the details of building mathematical optmization models** or the syntax of the Pyomo package. \n",
"3. To **automate many of the common steps of materials optimization,** speeding up the development of new models. \n",
"\n",
"As an example system, we will consider the minimization of cohesive energy in nanoclusters, recently demonstrated in:\n",
"\n",
"Isenberg, N. M., et al., \"Identification of Optimally Stable Nanocluster Geometries via Mathematical Optimization and Density Functional Theory,\" *Molecular Systems Design & Engineering* 5 (2020): 232-244. DOI: [10.1039/C9ME00108E](https://pubs.rsc.org/en/content/articlelanding/2020/me/c9me00108e#!divAbstract).\n",
"\n",
"We seek to identify the geometry that minimizes the cohesive energy of a nanocluster on the face-centered cubic (FCC) lattice. \n",
"As a model for cohesive energy, we use model based on the square-root of coordination number, refered to as the Tomanek model [[1]](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.28.665).\n",
"In the equation below, we define the normalized cohesive energy, as the normalized contribution of the square root of the coordination number. \n",
"\n",
"$$\\hat{E}^{\\text{surf}} = \\frac{1}{N \\sqrt{12}} \\displaystyle \\sum_i \\sqrt{CN_i} $$\n",
"\n",
"In the following sections, we demonstrate the basic approach for importing the MatOpt package, specifying the design space, formulating an optimization model, solving the optimization problem, and then outputting results. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importing Packages\n",
"\n",
"We start by importing several standard Python modules for convenience. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from math import sqrt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we import MatOpt."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from idaes.apps.matopt import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Representing Materials\n",
"To begin formulating a material optimization problem, we need several pieces of information about the design space.\n",
"Our goal is to generate a data structure for representing the choices in the design space, namely the choice of where to place atoms on FCC lattice sites. \n",
"\n",
"First, we define an ***FCCLattice*** object that holds the information about what sites should be included and which sites should be considered neighbors. \n",
"As argument to the lattice object, we are required to provide the interatomic distance. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"Lat = FCCLattice(IAD=2.770)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we define a ***Canvas*** object that links Cartesian coordinates to more abstract graph consisting of sites and neighbors.\n",
"We incrimentally construct a Canvas by first introducing a site at the origin of the coordinate system. \n",
"Then, we add \"two shells\" of neighbors, meaning that we introduce a shell of sites neighboring to the origin (12 for the FCC lattice) and then introduce another shell of neighbors to that group (42 additional sites, for a total of 55 sites). \n",
"The lattice object provides a *getNeighbors* method to identify these neighbors."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"Canv = Canvas()\n",
"Canv.addLocation(np.array([0,0,0],dtype=float))\n",
"Canv.addShells(2,Lat.getNeighbors)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we define a list of ***Atom*** objects to represent the building blocks of our materials. \n",
"We then use a ***Design*** object to represent the conjunction of a Canvas with a specific arrangement of building blocks.\n",
"The Design object will be used to represent the material decisions made during the solution of material optimization models. \n",
"\n",
"Before applying optimization, we can use the Design object to plot the sites of the Canvas and ensure that we constructed the intended design space. \n",
"We include several parsers to basic crystal structure file formats such as [XYZ](https://openbabel.org/docs/dev/FileFormats/XYZ_cartesian_coordinates_format.html), [PDB](https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/introduction), [POSCAR](https://www.vasp.at/wiki/index.php/POSCAR), and [CFG](http://li.mit.edu/Archive/Graphics/A/index.html#standard_CFG). "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"Atoms = [Atom('Pt')]\n",
"D = Design(Canv,Atom('Pt'))\n",
"D.toPDB('canvas_sites.pdb')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Building a Model\n",
"\n",
"To hold the materials information, we create a ***MatOptModel*** object. \n",
"This will hold information about the relevant Canvas, Atoms, and material conformations that may be present in a system. \n",
"Additionally, we define a parameter for the desired size of the cluster which will be utilized later by several methods. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"N = 22\n",
"m = MatOptModel(Canv,Atoms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The MatOptModel additionally hold lists of ***MaterialDescriptor*** objects that define the relevant material desriptors. \n",
"By default, several universal site descriptors are pre-defined in the model. \n",
"From these, all other material descriptors can be defined.\n",
"\n",
"| Descriptor | Explanation |\n",
"|:--------------|:------------|\n",
"| ***m.Yik*** | Presence of a building block of type k at site i |\n",
"| ***m.Yi*** | Presence of any type of building block at site i |\n",
"| ***m.Xijkl*** | Presence of a building block of type k at site i and a building block of type l at site j |\n",
"| ***m.Xij*** | Presence of any building block at site i and any building block at site j |\n",
"| ***m.Cikl*** | Count of neighbors of type l next to a building block of type k at site i |\n",
"| ***m.Ci*** | Count of any type of neighbors next to a building block at site i |\n",
"\n",
"User-specified descriptors are defined by ***DescriptorRule*** objects in conjunction with ***Expr*** expression objects. \n",
"Available expressions include:\n",
"\n",
"| Expression | Explanation |\n",
"|:-------------------------|:------------|\n",
"| ***LinearExpr*** | Multiplication and addition of coefficients to distinct MaterialDescriptors |\n",
"| ***SiteCombination*** | Summation of site contributions from two sites |\n",
"| ***SumNeighborSites*** | Summation of site contributions from all neighboring sites |\n",
"| ***SumNeighborBonds*** | Summation of bond contributions to all neighboring sites |\n",
"| ***SumSites*** | Summation across sites |\n",
"| ***SumBonds*** | Summation across bonds |\n",
"| ***SumSiteTypes*** | Summation across site types |\n",
"| ***SumBondTypes*** | Summation across bond types |\n",
"| ***SumSitesAndTypes*** | Summation across sites and site types |\n",
"| ***SumBondsAndTypes*** | Summation across bonds and bond types |\n",
"| ***SumConfs*** | Summation across conformation types |\n",
"| ***SumSitesAndConfs*** | Summation across sites and conformation types |\n",
"\n",
"Several types of DescriptorRules are available. \n",
"\n",
"| Rule | Explanation |\n",
"|:------------------------------|:------------|\n",
"| ***LessThan*** | Descriptor less than or equal to an expression |\n",
"| ***EqualTo*** | Descriptor equal to an expression |\n",
"| ***GreaterThan*** | Descriptor greater than or equal to an expression |\n",
"| ***FixedTo*** | Descriptor fixed to a scalar value |\n",
"| ***PiecewiseLinear*** | Descriptor equal to the evaluation of a piecewise linear function |\n",
"| ***Implies*** | Indicator descriptor that imposes other constraints if equal to 1 |\n",
"| ***NegImplies*** | Indicator descriptor that imposes other constraints if equal to 0 |\n",
"| ***ImpliesSiteCombination*** | Indicator bond-indexed descriptor that imposes constraints on the two sites |\n",
"| ***ImpliesNeighbors*** | Indicator site-indexed descriptor that imposes constraints on neighboring sites |\n",
"\n",
"From the combination of pre-defined descriptors, expressions, and rules we can specify a wide variety of other descriptors. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the context of nanocluster cohesive energy minimization, we would first like to define the square root of the coordination number. \n",
"We achieve this by calling the ***addSitesDescriptor*** method on MatOptModel, passing the information necessary to create a ***PiecewiseLinear*** rule to correctly define the square root of coordination at the integer coordination number values. \n",
"Note that we use ***m.Ci***, the pre-defined basic variable for the count of neighboring building blocks and equivalent to the coordination number in this system, as the argument for the piecewise linear function. \n",
"We use basic Python lists to express the data for the piecewise linear function values at integer numbers of coordination. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"Vals = [sqrt(CN) for CN in range(0,13)]\n",
"BPs = [CN for CN in range(0,13)]\n",
"m.addSitesDescriptor('CNRi',bounds=(0,sqrt(12)),integer=False,\n",
" rules=PiecewiseLinear(values=Vals,\n",
" breakpoints=BPs,\n",
" input_desc=m.Ci))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we define a global (i.e., not indexed by sites or bonds) descriptor for the cohesive energy of the nanocluster. \n",
"We us a simple ***EqualTo*** rule to set this descriptor equal to a normalized sum of contributions from the square root coordination number descriptor. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"m.addGlobalDescriptor('Ecoh',rules=EqualTo(SumSites(desc=m.CNRi,\n",
" coefs=(1/(N*sqrt(12))))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we create a descriptor for the size of the nanocluster. \n",
"We set bounds on this descriptor to effectively constrain the design space to only include clusters of the desired size, *N*. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"m.addGlobalDescriptor('Size',bounds=(N,N),\n",
" rules=EqualTo(SumSites(desc=m.Yi)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving the Model\n",
"Once all the relevant information in the model is provided, we search for optimal designs that maximize one of the descriptors. \n",
"In this example, we provide the descriptor for coehisver energy as the target functionality. \n",
"Additionally, we specify a time limit in seconds as a keyword argument to the maximize method. \n",
"For more information, see the documentation of the maximize function, available in the source code or by using the Python *help* function."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Help on function maximize in module matopt.opt.mat_modeling:\n",
"\n",
"maximize(self, func, **kwargs)\n",
" Method to maximize a target functionality of the material model.\n",
" \n",
" Args:\n",
" func (``MaterialDescriptor``/``Expr``): Material functionality to optimize.\n",
" **kwargs: Arguments to ``MatOptModel.optimize``\n",
" \n",
" Returns:\n",
" (``Design``/list<``Design``>) Optimal designs.\n",
" \n",
" Raises:\n",
" ``pyutilib.ApplicationError`` if MatOpt can not find usable solver (CPLEX or NEOS-CPLEX)\n",
" \n",
" See ``MatOptModel.optimize`` method for details.\n",
"\n",
"Help on function optimize in module matopt.opt.mat_modeling:\n",
"\n",
"optimize(self, func, sense, nSolns=1, tee=True, disp=1, keepfiles=False, tilim=3600, trelim=None, solver='cplex')\n",
" Method to create and optimize the materials design problem.\n",
" \n",
" This method automatically creates a new optimization model every \n",
" time it is called. Then, it solves the model via Pyomo with the \n",
" CPLEX solver.\n",
" \n",
" If multiple solutions (called a 'solution pool') are desired, then\n",
" the nSolns argument can be provided and the populate method will \n",
" be called instead. \n",
" \n",
" Args:\n",
" func (``MaterialDescriptor``/``Expr``): Material functionality to optimize.\n",
" sense (int): flag to indicate the choice to minimize or maximize the functionality of interest.\n",
" Choices: minimize/maximize (Pyomo constants 1,-1 respectively)\n",
" nSolns (int): Optional, number of Design objects to return.\n",
" Default: 1 (See ``MatOptModel.populate`` for more information)\n",
" tee (bool): Optional, flag to turn on solver output.\n",
" Default: True\n",
" disp (int): Optional, flag to control level of MatOpt output.\n",
" Choices: 0: No MatOpt output (other than solver tee) 1: MatOpt output for outer level method 2: MatOpt output for solution pool & individual solns.\n",
" Default: 1\n",
" keepfiles (bool): Optional, flag to save temporary pyomo files.\n",
" Default: True\n",
" tilim (float): Optional, solver time limit (in seconds).\n",
" Default: 3600\n",
" trelim (float): Optional, solver tree memeory limit (in MB).\n",
" Default: None (i.e., Pyomo/CPLEX default)\n",
" solver (str): Solver choice. Currently only cplex or neos-cplex are supported\n",
" Default: cplex\n",
" \n",
" Returns:\n",
" (``Design``/list<``Design``>) Optimal design or designs, depending on the number of solutions requested by argument ``nSolns``.\n",
" \n",
" Raises:\n",
" ``pyutilib.ApplicationError`` if MatOpt can not find usable solver (CPLEX or NEOS-CPLEX)\n",
"\n"
]
}
],
"source": [
"help(MatOptModel.maximize)\n",
"help(MatOptModel.optimize)"
]
},
{
"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 Community Edition 12.9.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\\tmp3lxwcq_l.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: 100\n",
"CPLEX> Problem 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmpdzztaw68.pyomo.lp' read.\n",
"Read time = 0.00 sec. (0.06 ticks)\n",
"CPLEX> Problem name : C:\\Users\\xiang\\AppData\\Local\\Temp\\tmpdzztaw68.pyomo.lp\n",
"Objective sense : Maximize\n",
"Variables : 426 [Nneg: 1, Fix: 1, Box: 13, Free: 157,\n",
" Binary: 241, General Integer: 13]\n",
"Objective nonzeros : 1\n",
"Linear constraints : 583 [Less: 515, Greater: 13, Equal: 55]\n",
" Nonzeros : 1866\n",
" RHS nonzeros : 86\n",
"\n",
"Variables : Min LB: 0.0000000 Max UB: 12.00000 \n",
"Objective nonzeros : Min : 1.000000 Max : 1.000000 \n",
"Linear constraints :\n",
" Nonzeros : Min : 0.04811252 Max : 12.00000 \n",
" RHS nonzeros : Min : 1.000000 Max : 1.000000 \n",
"CPLEX> CPXPARAM_TimeLimit 100\n",
"CPXPARAM_MIP_Tolerances_AbsMIPGap 0\n",
"CPXPARAM_MIP_Tolerances_MIPGap 0\n",
"Tried aggregator 2 times.\n",
"MIP Presolve eliminated 171 rows and 147 columns.\n",
"MIP Presolve modified 12 coefficients.\n",
"Aggregator did 38 substitutions.\n",
"Reduced MIP has 374 rows, 241 columns, and 1149 nonzeros.\n",
"Reduced MIP has 156 binaries, 0 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.02 sec. (1.49 ticks)\n",
"Found incumbent of value 0.521849 after 0.02 sec. (2.09 ticks)\n",
"Probing time = 0.00 sec. (1.36 ticks)\n",
"Cover probing fixed 0 vars, tightened 1 bounds.\n",
"Tried aggregator 1 time.\n",
"MIP Presolve eliminated 2 rows and 1 columns.\n",
"Reduced MIP has 372 rows, 240 columns, and 1144 nonzeros.\n",
"Reduced MIP has 156 binaries, 0 generals, 0 SOSs, and 0 indicators.\n",
"Presolve time = 0.00 sec. (1.50 ticks)\n",
"Probing time = 0.00 sec. (1.31 ticks)\n",
"Clique table members: 880.\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.00 sec. (3.34 ticks)\n",
"\n",
" Nodes Cuts/\n",
" Node Left Objective IInf Best Integer Best Bound ItCnt Gap\n",
"\n",
"* 0+ 0 0.5218 2.1667 315.19%\n",
" 0 0 0.9861 103 0.5218 0.9861 291 88.96%\n",
" 0 0 0.6455 96 0.5218 Cuts: 67 417 23.69%\n",
"* 0+ 0 0.5485 0.6455 17.68%\n",
" 0 0 0.6393 99 0.5485 Cuts: 159 460 16.55%\n",
" 0 0 0.6390 98 0.5485 ZeroHalf: 19 465 16.50%\n",
" 0 0 0.6356 99 0.5485 Cuts: 10 470 15.87%\n",
" 0 0 0.6336 97 0.5485 Cuts: 9 485 15.51%\n",
" 0 0 0.6108 97 0.5485 Cuts: 31 533 11.37%\n",
" 0 0 0.6088 98 0.5485 ZeroHalf: 20 548 10.99%\n",
" 0 0 0.6072 98 0.5485 ZeroHalf: 5 559 10.71%\n",
" 0 0 0.6032 97 0.5485 ZeroHalf: 22 579 9.97%\n",
" 0 0 0.6029 98 0.5485 Cuts: 20 588 9.92%\n",
" 0 0 0.5970 93 0.5485 ZeroHalf: 26 610 8.83%\n",
" 0 0 0.5970 96 0.5485 Cuts: 10 617 8.83%\n",
" 0 0 0.5969 95 0.5485 ZeroHalf: 7 627 8.82%\n",
"* 0+ 0 0.5500 0.5969 8.52%\n",
" 0 0 cutoff 0.5500 0.5500 627 0.00%\n",
"Elapsed time = 0.14 sec. (69.13 ticks, tree = 0.01 MB, solutions = 3)\n",
"\n",
"Clique cuts applied: 24\n",
"Implied bound cuts applied: 20\n",
"Zero-half cuts applied: 8\n",
"Lift and project cuts applied: 1\n",
"\n",
"Root node processing (before b&c):\n",
" Real time = 0.16 sec. (69.15 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.16 sec. (69.15 ticks)\n",
"\n",
"Solution pool: 3 solutions saved.\n",
"\n",
"MIP - Integer optimal solution: Objective = 5.5003296046e-01\n",
"Solution time = 0.16 sec. Iterations = 627 Nodes = 0\n",
"Deterministic time = 69.15 ticks (443.30 ticks/sec)\n",
"\n",
"CPLEX> Incumbent solution written to file 'C:\\Users\\xiang\\AppData\\Local\\Temp\\tmp72k3_0xk.cplex.sol'.\n",
"CPLEX> The solver exited normally.\n",
"A feasible and provably optimal solution is available.\n",
"The Design has objective: 0.5500329604578591\n"
]
}
],
"source": [
"D = None\n",
"try:\n",
" D = m.maximize(m.Ecoh,tilim=100)\n",
"except:\n",
" print('MaOpt can not find usable solver (CPLEX or NEOS-CPLEX)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Processing Results\n",
"\n",
"If a result is found, we can write it to file and plot with visualization software. \n",
"We provide interfaces to several standard crystal structure file formats, including [XYZ](https://openbabel.org/docs/dev/FileFormats/XYZ_cartesian_coordinates_format.html), [PDB](https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/introduction), [POSCAR](https://www.vasp.at/wiki/index.php/POSCAR), and [CFG](http://li.mit.edu/Archive/Graphics/A/index.html#standard_CFG). "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"if(D is not None):\n",
" D.toPDB('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.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}