The Kitchin Research Group

Almost 10 years ago I wrote about finding the minimum distance from a point to a curve using a constrained optimization. At that time, the way to do this used scipy.optimize.fmin_coblya. I learned today from a student, that sometimes this method fails! I reproduce the code here, updated for Python 3, some style updates, and to show it does indeed fail sometimes, notably when the point is "outside" the parabola.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla

def f(x):
    return x**2

for P in np.array([[0.5, 2],
                   [2, 2],
                   [-1, 2],
                   [-2, 2],
                   [0, 0.5],
                   [0, -0.5]]):
    
    def objective(X):
        X = np.array(X)
        return np.linalg.norm(X - P)

    def c1(X):
        x,y = X
        return f(x) - y

    X = fmin_cobyla(objective, x0=[P[0], f(P[0])], cons=[c1])

    print(f'The minimum distance is {objective(X):1.2f}. Constraint satisfied = {c1(X) < 1e-6}')

    # Verify the vector to this point is normal to the tangent of the curve
    # position vector from curve to point
    v1 = np.array(P) - np.array(X)
    # position vector
    v2 = np.array([1, 2.0 * X[0]])
    print('dot(v1, v2) = ', np.dot(v1, v2))

    x = np.linspace(-2, 2, 100)

    plt.plot(x, f(x), 'r-', label='f(x)')
    plt.plot(P[0], P[1], 'bo', label='point')
    plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')
    plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
    plt.axis('equal')
    plt.xlabel('x')
    plt.ylabel('y')    

The minimum distance is 0.86. Constraint satisfied = True dot(v1, v2) = 0.0002913487659186309 The minimum distance is 0.00. Constraint satisfied = False dot(v1, v2) = 0.00021460906432962284 The minimum distance is 0.39. Constraint satisfied = True dot(v1, v2) = 0.00014271520451364372 The minimum distance is 0.00. Constraint satisfied = False dot(v1, v2) = -0.0004089466778209598 The minimum distance is 0.50. Constraint satisfied = True dot(v1, v2) = 1.9999998429305957e-12 The minimum distance is 0.00. Constraint satisfied = False dot(v1, v2) = 8.588744170160093e-06

So, sure enough, the optimizer is failing to find a solution that meets the constraint. It is strange it does not work on the outside. That is almost certainly an algorithm problem. Here we solve it nearly identically with the more modern scipy.optimize.minimize function, and it converges every time.

from scipy.optimize import minimize

for P in np.array([[0.5, 2],
                   [2, 2],
                   [-1, 2],
                   [-2, 2],
                   [0, 0.5],
                   [0, -0.5]]):
    
    def objective(X):
        X = np.array(X)
        return np.linalg.norm(X - P)

    def c1(X):
        x,y = X
        return f(x) - y

    sol = minimize(objective, x0=[P[0], f(P[0])], constraints={'type': 'eq', 'fun': c1})
    X = sol.x

    print(f'The minimum distance is {objective(X):1.2f}. Constraint satisfied = {sol.status < 1e-6}')

    # Verify the vector to this point is normal to the tangent of the curve
    # position vector from curve to point
    v1 = np.array(P) - np.array(X)
    # position vector
    v2 = np.array([1, 2.0 * X[0]])
    print('dot(v1, v2) = ', np.dot(v1, v2))

    x = np.linspace(-2, 2, 100)

    plt.plot(x, f(x), 'r-', label='f(x)')
    plt.plot(P[0], P[1], 'bo', label='point')
    plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')
    plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
    plt.axis('equal')
    plt.xlabel('x')
    plt.ylabel('y')

The minimum distance is 0.86. Constraint satisfied = True dot(v1, v2) = 1.0701251773603815e-08 The minimum distance is 0.55. Constraint satisfied = True dot(v1, v2) = -0.0005793028003104883 The minimum distance is 0.39. Constraint satisfied = True dot(v1, v2) = -1.869272921939391e-05 The minimum distance is 0.55. Constraint satisfied = True dot(v1, v2) = 0.0005792953298950909 The minimum distance is 0.50. Constraint satisfied = True dot(v1, v2) = 0.0 The minimum distance is 0.50. Constraint satisfied = True dot(v1, v2) = 0.0

There is no wisdom in fixing the first problem, here I just tried a newer optimization method. Out of the box with default settings it just worked. I did learn the answer is sensitive to the initial guess, so it could make sense to sample the function and find the point that is closest as the initial guess, but here the simple heuristic guess I used worked fine.

Copyright (C) 2023 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.5.5

Check out the video at:

Nobody likes to run expensive jobs more than necessary, so cache solutions are often used where you save the results, and look them up later. There is functools.cache in Python, but it is memory only, and not persistent, so you start over in a new session.

For persistent cache, joblib (https://joblib.readthedocs.io/en/latest/) is a standard tool for this. Here is a simple example:

from joblib import Memory
location = '/Users/jkitchin/Dropbox/emacs/journal/2023/02/01/joblib_cache/joblib_cache'
memory = Memory(location, verbose=0)

@memory.cache
def fun(x=1.0):
    print('If you see this, go get a coffee while it runs 🐌.')
    return x**2

print(fun(2)) # Runs the function
print(fun(2)) # Looks up the cached value

If you see this, go get a coffee while it runs 🐌. 4 4

That works because joblib saves the results in a file in the location you specify.

Here is another example with another arg.

@memory.cache
def f2(x=1.0, a=3):
    print(f'If you see this, go get a coffee while it runs. a={"🐌"*a}')
    return a*x**2

(f2(2),       # Runs function
 f2(2, a=3),  # calls cache
 f2(2, 4))    # Runs another function because a changed

If you see this, go get a coffee while it runs. a=🐌🐌🐌 If you see this, go get a coffee while it runs. a=🐌🐌🐌🐌

Here, we look up from the cache each time.

(f2(2), f2(2, a=3), f2(2, 4))

a = 3

def f3(x=1.0):
    return a*x**2

f3(2)

a=0
f3(2)

a = 3
@memory.cache
def f4(x=1.0):
    print('If you see this, go get a coffee while it runs 🐌.')
    return a*x**2

(f4(2), f4(2), f4(2))

a = 0
(f4(2), f4(2), f4(2))

from ase.build import bulk
from ase.calculators.emt import EMT

atoms = bulk('Pd')
atoms.set_calculator(EMT())

@memory.cache
def f(atoms):
    print('If you see this, go get a coffee.')
    return atoms.get_potential_energy()

(f(atoms), f(atoms))

import joblib
atoms = bulk('Pd')
atoms.set_calculator(EMT())
joblib.hash(atoms), atoms._calc.results

atoms.get_potential_energy()
joblib.hash(atoms), atoms._calc.results

from joblib import Memory
location = '/Users/jkitchin/Dropbox/emacs/journal/2023/02/01/joblib_cache/joblib_cache'
memory = Memory(location, verbose=0)

a = 3
@memory.cache
def f4(x=1.0):
    print('If you see this, go get a coffee while it runs')
    return a*x**2

print((f4(2), f4(2), f4(2)))

@memory.cache
def f4(x=1.0):
    'add a ds.'
    # comment
    print('If you see this, go get a coffee while it runs')
    return a*x**2

print((f4(2), f4(2), f4(2)))

import pycse
pycse.__version__, pycse.__file__

from pycse.hashcache import hashcache
  
hashcache.location = "/Users/jkitchin/Dropbox/emacs/journal/2023/02/01/cache"
hashcache.verbose = False

@hashcache
def h1(x):
    print('This runs soo slow... Go get a coffee')
    return x**2

h1(2), h1(2)

a = 3
@hashcache
def h4(x=1.0):
    print('If you see this, go get a coffee while it runs')
    return a*x**2

(h4(2), h4(2), h4(2))

a=0
(h4(2), h4(2), h4(2))

from ase.build import bulk
from ase.calculators.emt import EMT

atoms = bulk('Pd')
atoms.set_calculator(EMT())

@hashcache
def h(atoms):
    print('If you see this, go get a coffee.')
    return atoms.get_potential_energy()

(h(atoms), h(atoms), h(atoms))

print(h1(2), h1(2))

from pycse.hashcache import hashcache
hashcache.location = "/Users/jkitchin/Dropbox/emacs/journal/2023/02/01/cache"

@hashcache
def h1(x):
    print('This runs soo slow... Go get a coffee')
    return x**2

print(h1(2), h1(2))

a = 3
@hashcache
def h4(x=1.0):
    'doc string'
    # comments
    print('If you see this, go get a coffee while it runs')
    return a*x**2    

(h4(2), h4(2), h4(2))

Copyright (C) 2023 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.5.5

2022 was an interesting one. I still did all of my teaching remotely, but spent more time going in to the office for meetings, and have gotten back to professional travel for meetings.

2. Publications

Our work this past year was divided in a few efforts. We had some work in method development, e.g. in uncertainty quantification, automatic differentiation, and vector search.

1. Yang, Y., Achar, S. K., & Kitchin, J. R. (2022). Evaluation of the degree of rate control via automatic differentiation. AIChE Journal, 68(6). http://dx.doi.org/10.1002/aic.17653
2. Yang, Y., Liu, M., & Kitchin, J. R. (2022). Neural network embeddings based similarity search method for atomistic systems. Digital Discovery. http://dx.doi.org/10.1039/d2dd00055e
3. Zhan, N., & Kitchin, J. R. (2022). Model-specific to model-general uncertainty for physical properties. Industrial & Engineering Chemistry Research, 1–04706. http://dx.doi.org/10.1021/acs.iecr.1c04706

This collaborative work on large catalyst models is was especially exciting. Stay tuned for many advances in this area in 2023.

1. Kolluru, A., Shuaibi, M., Palizhati, A., Shoghi, N., Das, A., Wood, B., Zitnick, C. L., … (2022). Open challenges in developing generalizable large-scale machine-learning models for catalyst discovery. ACS Catalysis, 12(14), 8572–8581. http://dx.doi.org/10.1021/acscatal.2c02291

We also wrote some collaborative papers on our work in high-throughput discovery of hydrogen evolution catalysts and segregation.

1. Broderick, K., Lopato, E., Wander, B., Bernhard, S., Kitchin, J., & Ulissi, Z. (2022). Identifying limitations in screening high-throughput photocatalytic bimetallic nanoparticles with machine-learned hydrogen adsorptions. Applied Catalysis B: Environmental, 121959. http://dx.doi.org/10.1016/j.apcatb.2022.121959
2. Bhat, M., Lopato, E., Simon, Z. C., Millstone, J. E., Bernhard, S., & Kitchin, J. R. (2022). Accelerated optimization of pure metal and ligand compositions for light-driven hydrogen production. Reaction Chemistry & Engineering. http://dx.doi.org/10.1039/d1re00441g
3. Yilin Yang, Zhitao Guo, Andrew Gellman and John Kitchin, Simulating segregation in a ternary Cu-Pd-Au alloy with density functional theory, machine learning and Monte Carlo simulations, J. Phys. Chem. C, 126, 4, 1800-1808. (2022). https://pubs.acs.org/doi/abs/10.1021/acs.jpcc.1c09647

It was another big year in citations for us (https://scholar.google.com/citations?user=jD_4h7sAAAAJ)!

Copyright (C) 2023 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.5.5

Hydrogen adsorption energies have long been used in screening to identify promising hydrogen evolution catalysts. In this work we combine a high-throughput experimental study of 5300 different bimetallic catalysts with a high-throughput computational screen of 16M adsorption energies to see how well adsorption energies work for this purpose. We developed a workflow to combine these data sets that accounts for surface stability and adsorption site distributions on alloy surfaces. We find that thermodynamically favorable adsorption energies are necessary to observe high activity, but they are not sufficient, and do not always lead to high activity.

@article{broderick-2022-ident-limit,
  author = {Kirby Broderick and Eric Lopato and Brook Wander and Stefan Bernhard and John Kitchin and Zachary Ulissi},
  title = {Identifying Limitations in Screening High-Throughput Photocatalytic Bimetallic Nanoparticles With Machine-Learned Hydrogen Adsorptions},
  journal = {Applied Catalysis B: Environmental},
  volume = {nil},
  number = {nil},
  pages = {121959},
  year = {2022},
  doi = {10.1016/j.apcatb.2022.121959},
  url = {http://dx.doi.org/10.1016/j.apcatb.2022.121959},
  DATE_ADDED = {Thu Sep 22 07:46:56 2022},
}

Copyright (C) 2022 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.5.1

Searching for atomic structures in databases is like finding a needle in the haystack. It is difficult to construct a query that finds what you want, without finding nothing, or everything! It is difficult to use atomic coordinates because they are sensitive to translations, rotations and permutations. There are many ways to construct equivalent unit cells that also make it difficult to uniquely query materials.

In this paper we show how to construct queries for atomic structures that allow you to quickly find similar atomic structures. We achieve this by using invariant fingerprint vectors from machine learning models coupled with approximate nearest neighbor vector search algorithms. We apply it to molecules, bulk materials and adsorbates on surfaces. We show how the geometric similarity in found atomic systems leads to better data sets for building new machine learning models, and that the found systems tend to show geometric and electronic structure similarity.

You can read more about this work here (it is Open Access): Yilin Yang, Mingie Liu, John Kitchin, Digital Discovery, 2022, https://doi.org/10.1039/D2DD00055E

@article{yang-2022-neural-networ,
  author =       {Yilin Yang and Mingjie Liu and John R. Kitchin},
  title =        {Neural Network Embeddings Based Similarity Search Method for
                  Atomistic Systems},
  journal =      {Digital Discovery},
  volume =       {},
  number =       {},
  pages =        {},
  year =         2022,
  doi =          {10.1039/d2dd00055e},
  url =          {https://doi.org/10.1039/D2DD00055E},
  DATE_ADDED =   {Mon Sep 12 17:21:30 2022},
}

Copyright (C) 2022 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.5.1