Bimetallic Nanocluster Design¶
Bimetallic Nanocluster Cohesive Energy Minimization - via Labeling¶
This notebook serves as an example application of the MatOpt framework. We consider an example optimization problem of identifying the global energy minimum bimetallic nanocluster configuration.
This is a continuation of the example given in *Monometallic_Nanocluster_Design.ipynb*. In this example, we will show how a very similar model can be used to optimize a bimetallic cluster by "labelling" the sites of a pre-defined monometallic cluster.
For more information, see: Yin, X., et al. "Designing Stable Bimetallic Nanoclusters via an Iterative Two-Step Optimization Approach." Molecular Systems Design & Engineering (2021), Accepted for Publication.
Importing Packages¶
We start by importing several standard Python modules for convienience.
import numpy as np
from math import sqrt
Then, we import the MatOpt package in its entirety.
from idaes.apps.matopt import *
Setting Up a Material System¶
We first identify the optimal metal-independent nanocluster shape, using the code that wsas demonstrated in Monometallic_Nanocluster_Design.ipynb.
Lat = FCCLattice(IAD=1.0)
Canv = Canvas()
Canv.addLocation(np.array([0,0,0],dtype=float))
Canv.addShells(1,Lat.getNeighbors)
Atoms = [Atom('Cu')]
N = 6
m = MatOptModel(Canv,Atoms)
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))
m.addGlobalDescriptor('Ecoh',rules=EqualTo(SumSites(desc=m.CNRi,
coefs=(1/(N*sqrt(12))))))
m.addGlobalDescriptor('Size',bounds=(N,N),
rules=EqualTo(SumSites(desc=m.Yi)))
# This step requires the CPLEX solver
try:
D = m.maximize(m.Ecoh,tilim=100)
except Exception as err:
D = None # rest of this notebook won't do much
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/tmpt6868vkb.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/tmp5e789y_h.pyomo.lp' read.
Read time = 0.00 sec. (0.06 ticks)
CPLEX> Problem name : /tmp/tmp5e789y_h.pyomo.lp
Objective sense : Maximize
Variables : 426 [Nneg: 1, Fix: 1, Box: 13, Free: 157,
Binary: 241, General Integer: 13]
Objective nonzeros : 1
Linear constraints : 583 [Less: 515, Greater: 13, Equal: 55]
Nonzeros : 1866
RHS nonzeros : 86
Variables : Min LB: 0.000000 Max UB: 12.00000
Objective nonzeros : Min : 1.000000 Max : 1.000000
Linear constraints :
Nonzeros : Min : 0.04811252 Max : 12.00000
RHS nonzeros : Min : 1.000000 Max : 1.000000
CPLEX> CPXPARAM_TimeLimit 100
CPXPARAM_MIP_Tolerances_AbsMIPGap 0
CPXPARAM_MIP_Tolerances_MIPGap 0
Tried aggregator 2 times.
MIP Presolve eliminated 171 rows and 147 columns.
MIP Presolve modified 12 coefficients.
Aggregator did 38 substitutions.
Reduced MIP has 374 rows, 241 columns, and 1149 nonzeros.
Reduced MIP has 156 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (1.49 ticks)
Found incumbent of value 0.521849 after 0.00 sec. (2.09 ticks)
Probing time = 0.00 sec. (1.36 ticks)
Cover probing fixed 0 vars, tightened 1 bounds.
Tried aggregator 1 time.
MIP Presolve eliminated 2 rows and 1 columns.
Reduced MIP has 372 rows, 240 columns, and 1144 nonzeros.
Reduced MIP has 156 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (1.50 ticks)
Probing time = 0.00 sec. (1.31 ticks)
Clique table members: 880.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 32 threads.
Root relaxation solution time = 0.00 sec. (3.19 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
* 0+ 0 0.5218 2.1667 315.19%
0 0 0.9861 105 0.5218 0.9861 278 88.96%
0 0 0.6455 96 0.5218 Cuts: 68 386 23.69%
* 0+ 0 0.5485 0.6455 17.68%
0 0 0.6351 99 0.5485 Cuts: 82 430 15.80%
0 0 0.6329 93 0.5485 ZeroHalf: 24 446 15.39%
0 0 0.6326 99 0.5485 ZeroHalf: 34 465 15.33%
0 0 0.6324 99 0.5485 Cuts: 22 476 15.29%
0 0 0.6067 96 0.5485 Cuts: 41 533 10.62%
0 0 0.6050 96 0.5485 ZeroHalf: 19 553 10.30%
0 0 0.6048 97 0.5485 ZeroHalf: 28 563 10.27%
0 0 0.6033 99 0.5485 ZeroHalf: 18 591 9.99%
0 0 0.6021 97 0.5485 ZeroHalf: 20 619 9.77%
0 0 0.6005 100 0.5485 ZeroHalf: 41 637 9.48%
0 0 0.5994 95 0.5485 ZeroHalf: 35 671 9.28%
0 0 0.5988 98 0.5485 ZeroHalf: 12 686 9.17%
0 0 0.5983 99 0.5485 ZeroHalf: 19 719 9.08%
0 0 0.5980 99 0.5485 ZeroHalf: 13 731 9.02%
0 0 0.5976 99 0.5485 ZeroHalf: 11 753 8.95%
0 0 0.5973 94 0.5485 ZeroHalf: 15 769 8.90%
0 0 0.5968 99 0.5485 ZeroHalf: 4 793 8.81%
0 0 0.5963 99 0.5485 ZeroHalf: 4 815 8.72%
* 0+ 0 0.5500 0.5963 8.42%
0 0 cutoff 0.5500 0.5500 831 0.00%
Elapsed time = 0.21 sec. (98.42 ticks, tree = 0.01 MB, solutions = 3)
Clique cuts applied: 10
Implied bound cuts applied: 20
Zero-half cuts applied: 17
Lift and project cuts applied: 1
Gomory fractional cuts applied: 1
Root node processing (before b&c):
Real time = 0.21 sec. (98.44 ticks)
Parallel b&c, 32 threads:
Real time = 0.00 sec. (0.00 ticks)
Sync time (average) = 0.00 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.21 sec. (98.44 ticks)
Solution pool: 3 solutions saved.
MIP - Integer optimal solution: Objective = 5.5003296046e-01
Solution time = 0.22 sec. Iterations = 831 Nodes = 0
Deterministic time = 98.45 ticks (444.83 ticks/sec)
CPLEX> Incumbent solution written to file '/tmp/tmp06rzq3xj.cplex.sol'.
CPLEX> The solver exited normally.
A feasible and provably optimal solution is available.
The Design has objective: 0.5500329604578591
We take the locations from the optimal monometallic problem to initialize a *Canvas* object for the bimetallic case.
Canv = Canvas()
if D: # may be None if CPLEX was not found
for i in range(len(D)):
if(D.Contents[i] is not None):
Canv.addLocation(D.Canvas.Points[i])
Canv.setNeighborsFromFunc(Lat.getNeighbors)
Additionally, we create a few data structures for holding bimetallic material information. First, we make a list of multiple *Atom* objects that will be the building blocks of the model. Next, we specify a dictionary with the bounds to impose on composition.
Atoms = [Atom('Cu'),Atom('Ag')]
CompBounds = {Atom('Cu'):(3,3),
Atom('Ag'):(3,3)}
Specifying an Optimization Model¶
We start by creating a *MatOptModel* object that will hold the information about the problem variables and constraints. At a minimum, ever model requires a Canvas object to be defined. Additionally, the list of building blocks and conformations that are present in the model should be defined.
m = MatOptModel(Canv,Atoms)
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.
To start, we inidcate that the choice to place an atom is fixed so that each canvas site is required to have an atom. This simplifies the problem significantly and results in a model that will seek to find the optimal labeling of metals on the nanocluster.
m.Yi.rules.append(FixedTo(1.0))
Next, we define a descriptor for the energy of bonds as a function of properties at each site. Since the locations of the atoms are fixed, the only decision is how to label each site as either Atom A or Atom B. This allows us to simplify the model and compute coefficients that rely on coordination number. In the block below, we implement the bimetallic model for bond energy defined in Yan et al., 2018.
GklCoefs = {(Atom('Cu'),Atom('Cu')):3.520,
(Atom('Cu'),Atom('Ag')):2.112,
(Atom('Ag'),Atom('Ag')):2.580,
(Atom('Ag'),Atom('Cu')):3.612}
BEijCoefs = {}
for i in range(len(Canv)):
CNi = sum(1 for _ in Canv.NeighborhoodIndexes[i] if _ is not None)
for j in Canv.NeighborhoodIndexes[i]:
if(j is not None):
CNj = sum(1 for _ in Canv.NeighborhoodIndexes[j] if _ is not None)
for k in Atoms:
for l in Atoms:
BEijCoefs[i,j,k,l] = GklCoefs[k,l]*1/sqrt(CNi) + GklCoefs[l,k]*1/sqrt(CNj)
m.addBondsDescriptor('BEij',
rules=EqualTo(SumBondTypes(m.Xijkl,coefs=BEijCoefs)),
symmetric_bonds=True)
Next, we define the cohesive energy as a sum of contributions from all BEij bond descriptors.
m.addGlobalDescriptor('Ecoh',rules=EqualTo(SumBonds(desc=m.BEij,
coefs=1/(N*sqrt(12)))))
Finally, we add constraints on the size and composition of the resulting designs.
m.addGlobalTypesDescriptor('Composition',bounds=CompBounds,
rules=EqualTo(SumSites(desc=m.Yik)))
Solving the Model¶
Once the model is fully specified, we can optimize in light of a global descriptor. In this example, we choose to maximize the cohesive energy defined previously. Additionally, we can specify basic optimization parameters such as the time limit and memory limit* for the optimizer.
D = None
try:
D = m.maximize(m.Ecoh,tilim=360,trelim=4096)
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/tmp_n7ozdru.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> New value for upper limit on size of tree in megabytes: 4096
CPLEX> Problem '/tmp/tmpdl2zsnb3.pyomo.lp' read.
Read time = 0.00 sec. (0.01 ticks)
CPLEX> Problem name : /tmp/tmpdl2zsnb3.pyomo.lp
Objective sense : Maximize
Variables : 71 [Nneg: 1, Fix: 2, Free: 12, Binary: 56]
Objective nonzeros : 1
Linear constraints : 159 [Less: 138, Equal: 21]
Nonzeros : 414
RHS nonzeros : 57
Variables : Min LB: 0.000000 Max UB: 3.000000
Objective nonzeros : Min : 1.000000 Max : 1.000000
Linear constraints :
Nonzeros : Min : 0.04811252 Max : 4.064546
RHS nonzeros : Min : 1.000000 Max : 1.000000
CPLEX> CPXPARAM_TimeLimit 360
CPXPARAM_MIP_Tolerances_AbsMIPGap 0
CPXPARAM_MIP_Tolerances_MIPGap 0
CPXPARAM_MIP_Limits_TreeMemory 4096
Tried aggregator 2 times.
MIP Presolve eliminated 20 rows and 15 columns.
Aggregator did 6 substitutions.
Reduced MIP has 133 rows, 50 columns, and 314 nonzeros.
Reduced MIP has 50 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.44 ticks)
Found incumbent of value 1.546227 after 0.00 sec. (0.80 ticks)
Probing time = 0.00 sec. (0.23 ticks)
Tried aggregator 1 time.
Reduced MIP has 133 rows, 50 columns, and 314 nonzeros.
Reduced MIP has 50 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.30 ticks)
Probing time = 0.00 sec. (0.23 ticks)
Clique table members: 263.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 32 threads.
Root relaxation solution time = 0.00 sec. (0.30 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
* 0+ 0 1.5462 6.5036 320.61%
* 0+ 0 1.6517 6.5036 293.74%
0 0 3.2518 50 1.6517 3.2518 50 96.87%
* 0+ 0 1.6556 3.2518 96.41%
0 0 2.1027 38 1.6556 Cuts: 112 96 27.01%
* 0+ 0 1.6754 2.1027 25.50%
0 0 1.8260 31 1.6754 Cuts: 41 120 8.99%
0 0 1.8217 16 1.6754 ZeroHalf: 9 121 8.73%
0 0 1.8121 33 1.6754 ZeroHalf: 2 124 8.16%
0 0 1.7548 15 1.6754 Cuts: 7 126 4.74%
0 0 cutoff 1.6754 129 ---
Elapsed time = 0.08 sec. (8.81 ticks, tree = 0.01 MB, solutions = 4)
Clique cuts applied: 40
Implied bound cuts applied: 6
Zero-half cuts applied: 13
Lift and project cuts applied: 1
Root node processing (before b&c):
Real time = 0.08 sec. (8.81 ticks)
Parallel b&c, 32 threads:
Real time = 0.00 sec. (0.00 ticks)
Sync time (average) = 0.00 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.08 sec. (8.81 ticks)
Solution pool: 4 solutions saved.
MIP - Integer optimal solution: Objective = 1.6754118151e+00
Solution time = 0.08 sec. Iterations = 129 Nodes = 0
Deterministic time = 8.81 ticks (106.68 ticks/sec)
CPLEX> Incumbent solution written to file '/tmp/tmplypn1ef_.cplex.sol'.
CPLEX> The solver exited normally.
A feasible and provably optimal solution is available.
The Design has objective: 1.6754118151174977
Processing Solutions¶
If a design was identified (optimal or otherwise), then a *Design* object is returned from the optimization method. The optimal design can be plotted via any of the supported parsers.
if D is not None:
D.toPDB('result.pdb')