Monometallic Nanocluster Design

Jupyter Notebook source file

Monometallic Nanocluster Design¶

In this module, we introduce the MatOpt interface for representing material properties and specifying optimization problems.

We have designed the interface with severl goals in mind:

1. To simplify the representation of nanostructured materials, streamlining the creation of materials optimization problems.
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.
3. To automate many of the common steps of materials optimization, speeding up the development of new models.

As an example system, we will consider the minimization of cohesive energy in nanoclusters, recently demonstrated in:

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.

We seek to identify the geometry that minimizes the cohesive energy of a nanocluster on the face-centered cubic (FCC) lattice. As a model for cohesive energy, we use model based on the square-root of coordination number, refered to as the Tomanek model [1]. In the equation below, we define the normalized cohesive energy, as the normalized contribution of the square root of the coordination number.

$$\hat{E}^{\text{surf}} = \frac{1}{N \sqrt{12}} \displaystyle \sum_i \sqrt{CN_i} $$

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.

Importing Packages¶

We start by importing several standard Python modules for convenience.

In [1]:

import numpy as np
from math import sqrt

Next, we import MatOpt.

In [2]:

from idaes.apps.matopt import *

Representing Materials¶

To begin formulating a material optimization problem, we need several pieces of information about the design space. 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.

First, we define an *FCCLattice* object that holds the information about what sites should be included and which sites should be considered neighbors. As argument to the lattice object, we are required to provide the interatomic distance.

In [3]:

Lat = FCCLattice(IAD=2.770)

Next, we define a *Canvas* object that links Cartesian coordinates to more abstract graph consisting of sites and neighbors. We incrimentally construct a Canvas by first introducing a site at the origin of the coordinate system. 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). The lattice object provides a getNeighbors method to identify these neighbors.

In [4]:

Canv = Canvas()
Canv.addLocation(np.array([0,0,0],dtype=float))
Canv.addShells(2,Lat.getNeighbors)

Finally, we define a list of *Atom* objects to represent the building blocks of our materials. We then use a *Design* object to represent the conjunction of a Canvas with a specific arrangement of building blocks. The Design object will be used to represent the material decisions made during the solution of material optimization models.

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. We include several parsers to basic crystal structure file formats such as XYZ, PDB, POSCAR, and CFG.

In [5]:

Atoms = [Atom('Pt')]
D = Design(Canv,Atom('Pt'))
D.toPDB('canvas_sites.pdb')

Building a Model¶

To hold the materials information, we create a *MatOptModel* object. This will hold information about the relevant Canvas, Atoms, and material conformations that may be present in a system. Additionally, we define a parameter for the desired size of the cluster which will be utilized later by several methods.

In [6]:

N = 22
m = MatOptModel(Canv,Atoms)

The MatOptModel additionally hold lists of *MaterialDescriptor* objects that define the relevant material desriptors. By default, several universal site descriptors are pre-defined in the model. From these, all other material descriptors can be defined.

Descriptor	Explanation
*m.Yik*	Presence of a building block of type k at site i
*m.Yi*	Presence of any type of building block at site i
*m.Xijkl*	Presence of a building block of type k at site i and a building block of type l at site j
*m.Xij*	Presence of any building block at site i and any building block at site j
*m.Cikl*	Count of neighbors of type l next to a building block of type k at site i
*m.Ci*	Count of any type of neighbors next to a building block at site i

User-specified descriptors are defined by *DescriptorRule* objects in conjunction with *Expr* expression objects. Available expressions include:

Expression	Explanation
*LinearExpr*	Multiplication and addition of coefficients to distinct MaterialDescriptors
*SiteCombination*	Summation of site contributions from two sites
*SumNeighborSites*	Summation of site contributions from all neighboring sites
*SumNeighborBonds*	Summation of bond contributions to all neighboring sites
*SumSites*	Summation across sites
*SumBonds*	Summation across bonds
*SumSiteTypes*	Summation across site types
*SumBondTypes*	Summation across bond types
*SumSitesAndTypes*	Summation across sites and site types
*SumBondsAndTypes*	Summation across bonds and bond types
*SumConfs*	Summation across conformation types
*SumSitesAndConfs*	Summation across sites and conformation types

Several types of DescriptorRules are available.

Rule	Explanation
*LessThan*	Descriptor less than or equal to an expression
*EqualTo*	Descriptor equal to an expression
*GreaterThan*	Descriptor greater than or equal to an expression
*FixedTo*	Descriptor fixed to a scalar value
*PiecewiseLinear*	Descriptor equal to the evaluation of a piecewise linear function
*Implies*	Indicator descriptor that imposes other constraints if equal to 1
*NegImplies*	Indicator descriptor that imposes other constraints if equal to 0
*ImpliesSiteCombination*	Indicator bond-indexed descriptor that imposes constraints on the two sites
*ImpliesNeighbors*	Indicator site-indexed descriptor that imposes constraints on neighboring sites

From the combination of pre-defined descriptors, expressions, and rules we can specify a wide variety of other descriptors.

In the context of nanocluster cohesive energy minimization, we would first like to define the square root of the coordination number. 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. 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. We use basic Python lists to express the data for the piecewise linear function values at integer numbers of coordination.

In [7]:

Vals = [sqrt(CN) for CN in range(0,13)]
BPs = [CN for CN in range(0,13)]
m.addSitesDescriptor('CNRi',bounds=(0,sqrt(12)),integer=False,
                     rules=PiecewiseLinear(values=Vals,
                                           breakpoints=BPs,
                                           input_desc=m.Ci))

Next, we define a global (i.e., not indexed by sites or bonds) descriptor for the cohesive energy of the nanocluster. We us a simple *EqualTo* rule to set this descriptor equal to a normalized sum of contributions from the square root coordination number descriptor.

In [8]:

m.addGlobalDescriptor('Ecoh',rules=EqualTo(SumSites(desc=m.CNRi,
                                                    coefs=(1/(N*sqrt(12))))))

Finally, we create a descriptor for the size of the nanocluster. We set bounds on this descriptor to effectively constrain the design space to only include clusters of the desired size, N.

In [9]:

m.addGlobalDescriptor('Size',bounds=(N,N),
                      rules=EqualTo(SumSites(desc=m.Yi)))

Solving the Model¶

Once all the relevant information in the model is provided, we search for optimal designs that maximize one of the descriptors. In this example, we provide the descriptor for coehisver energy as the target functionality. Additionally, we specify a time limit in seconds as a keyword argument to the maximize method. For more information, see the documentation of the maximize function, available in the source code or by using the Python help function.

In [10]:

help(MatOptModel.maximize)
help(MatOptModel.optimize)

Help on function maximize in module matopt.opt.mat_modeling:

maximize(self, func, **kwargs)
    Method to maximize a target functionality of the material model.
    
    Args:
        func (``MaterialDescriptor``/``Expr``): Material functionality to optimize.
        **kwargs: Arguments to ``MatOptModel.optimize``
    
    Returns:
        (``Design``/list<``Design``>) Optimal designs.
    
    Raises:
        ``pyomo.common.errors.ApplicationError`` if MatOpt can not find
        usable solver (CPLEX or NEOS-CPLEX)
    
    See ``MatOptModel.optimize`` method for details.

Help on function optimize in module matopt.opt.mat_modeling:

optimize(self, func, sense, nSolns=1, tee=True, disp=1, keepfiles=False, tilim=3600, trelim=None, solver='cplex')
    Method to create and optimize the materials design problem.
    
    This method automatically creates a new optimization model every
    time it is called. Then, it solves the model via Pyomo with the
    CPLEX solver.
    
    If multiple solutions (called a 'solution pool') are desired, then
    the nSolns argument can be provided and the populate method will
    be called instead.
    
    Args:
        func (``MaterialDescriptor``/``Expr``): Material functionality to optimize.
        sense (int): flag to indicate the choice to minimize or maximize the
            functionality of interest.
            Choices: minimize/maximize (Pyomo constants 1,-1 respectively)
        nSolns (int): Optional, number of Design objects to return.
            Default: 1 (See ``MatOptModel.populate`` for more information)
        tee (bool): Optional, flag to turn on solver output.
            Default: True
        disp (int): Optional, flag to control level of MatOpt output.
            Choices: 0: No MatOpt output (other than solver tee) 1: MatOpt
            output for outer level method 2: MatOpt output for solution pool &
            individual solns. Default: 1
        keepfiles (bool): Optional, flag to save temporary pyomo files.
            Default: True
        tilim (float): Optional, solver time limit (in seconds).
            Default: 3600
        trelim (float): Optional, solver tree memeory limit (in MB).
            Default: None (i.e., Pyomo/CPLEX default)
        solver (str): Solver choice. Currently only cplex or neos-cplex are supported
            Default: cplex
    
    Returns:
        (``Design``/list<``Design``>) Optimal design or designs, depending
        on the number of solutions requested by argument ``nSolns``.
    
    Raises:
        ``pyomo.common.errors.ApplicationError`` if MatOpt can not find
        usable solver (CPLEX or NEOS-CPLEX)

In [11]:

D = None
try:
    D = m.maximize(m.Ecoh,tilim=100)
except:
    print('MaOpt can not find usable solver (CPLEX or NEOS-CPLEX)')

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/tmp0a2sguf3.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: 100
CPLEX> Problem '/tmp/tmp1rkwydoq.pyomo.lp' read.
Read time = 0.01 sec. (0.28 ticks)
CPLEX> Problem name         : /tmp/tmp1rkwydoq.pyomo.lp
Objective sense      : Maximize
Variables            :    1920  [Nneg: 1,  Fix: 1,  Box: 55,  Free: 661,
                                 Binary: 1147,  General Integer: 55]
Objective nonzeros   :       1
Linear constraints   :    2839  [Less: 2561,  Greater: 55,  Equal: 223]
  Nonzeros           :    8904
  RHS nonzeros       :     488

Variables            : Min LB: 0.000000         Max UB: 22.00000       
Objective nonzeros   : Min   : 1.000000         Max   : 1.000000       
Linear constraints   :
  Nonzeros           : Min   : 0.01312160       Max   : 12.00000       
  RHS nonzeros       : Min   : 1.000000         Max   : 1.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 Results¶

If a result is found, we can write it to file and plot with visualization software. We provide interfaces to several standard crystal structure file formats, including XYZ, PDB, POSCAR, and CFG.

In [12]:

if(D is not None):
    D.toPDB('result.pdb')