Predicting Reactivity via Condensed Fukui Functions¶
Use AIMNet2 to compare how electron-donating and electron-withdrawing substituents change the local reactivity of a vinyl group.
Predicting Reactivity via Condensed Fukui Functions¶
This example uses the AIMNet2-NSE machine learning potential to compare a series of molecules with the common vinyl fragment \(\mathrm{CH_2{=}CH{-}R}\). The substituent \(R\) ranges from electron donating to electron withdrawing. This will demonstrate how the terminal vinyl carbon changes from being predominantly nucleophilic to predominantly electrophilic in nature.
For an atom \(A\), the condensed Fukui functions describe its response to adding or removing an electron. \(f_A^+\) measures the response to electron addition and is large at electrophilic sites, which accept electron density from nucleophiles. \(f_A^-\) measures the response to electron removal and is large at nucleophilic sites, which donate electron density to electrophiles.
The dual descriptor is \(\Delta f_A=f_A^+-f_A^-\). Positive values indicate predominantly electrophilic character, while negative values indicate predominantly nucleophilic character.
AIMNet2 evaluates the electron-added and electron-removed states at the optimized neutral geometry. These same calculations provide the vertical electron affinity and ionization energy.
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.axes import Axes
from scm.base import ChemicalSystem, Units
from scm.plams import AMSJob, Settings, plot_image_grid, view, view_atomic_property
Build the substituted alkenes¶
The SMILES strings are written from the terminal vinyl carbon, making it atom 1 in every molecule. Methoxy is a resonance donor, methyl is a weak donor through hyperconjugation, and cyano and formyl are electron-withdrawing groups. Ethene provides the unsubstituted reference.
substituted_alkenes = {
"Methyl vinyl ether": {"smiles": "C=COC", "group": r"$-\mathrm{OCH_3}$"},
"Propene": {"smiles": "C=CC", "group": r"$-\mathrm{CH_3}$"},
"Ethene": {"smiles": "C=C", "group": r"$-\mathrm{H}$"},
"Acrylonitrile": {"smiles": "C=CC#N", "group": r"$-\mathrm{CN}$"},
"Acrolein": {"smiles": "C=CC=O", "group": r"$-\mathrm{CHO}$"},
}
molecules = {
name: ChemicalSystem.from_smiles(details["smiles"])
for name, details in substituted_alkenes.items()
}
plot_image_grid(
{name: view(molecule, direction="along_pca3") for name, molecule in molecules.items()}, rows=1
);
Optimize the neutral structures and calculate the descriptors¶
When AMS Properties%Other is requested, the condensed Fukui functions are automatically calculated for an aperiodic ML potential that accepts total charge and returns individual atomic charges. This automatic behavior can be controlled with Fukui%Enabled. During a geometry optimization, the charged states are evaluated only after optimizing the neutral system. Their geometry is held fixed, so all reported energy differences are vertical.
def reactivity_settings() -> Settings:
settings = Settings()
settings.input.ams.Task = "GeometryOptimization"
settings.input.MLPotential.Model = "AIMNet2-NSE"
settings.runscript.nproc = 1
return settings
jobs = {}
for name, molecule in molecules.items():
job_name = name.lower().replace(" ", "_")
jobs[name] = AMSJob(name=job_name, molecule=molecule, settings=reactivity_settings())
jobs[name].run();
[24.09|12:00:05] JOB methyl_vinyl_ether STARTED
[24.09|12:00:05] JOB methyl_vinyl_ether RUNNING
[24.09|12:00:15] JOB methyl_vinyl_ether FINISHED
[24.09|12:00:15] JOB methyl_vinyl_ether SUCCESSFUL
[24.09|12:00:15] JOB propene STARTED
[24.09|12:00:15] JOB propene RUNNING
[24.09|12:00:24] JOB propene FINISHED
[24.09|12:00:24] JOB propene SUCCESSFUL
[24.09|12:00:24] JOB ethene STARTED
[24.09|12:00:24] JOB ethene RUNNING
[24.09|12:00:33] JOB ethene FINISHED
[24.09|12:00:33] JOB ethene SUCCESSFUL
[24.09|12:00:33] JOB acrylonitrile STARTED
[24.09|12:00:33] JOB acrylonitrile RUNNING
[24.09|12:00:41] JOB acrylonitrile FINISHED
[24.09|12:00:41] JOB acrylonitrile SUCCESSFUL
[24.09|12:00:41] JOB acrolein STARTED
[24.09|12:00:41] JOB acrolein RUNNING
[24.09|12:00:50] JOB acrolein FINISHED
[24.09|12:00:50] JOB acrolein SUCCESSFUL
Read and plot descriptor results¶
The Fukui arrays are stored in the MLPotential engine results. The helper below accepts a zero-based Python atom index. Thus, index 0 corresponds to atom 1: the terminal \(\mathrm{CH_2}\) carbon in every structure.
def read_reactivity_descriptors(job: AMSJob, atom_index: int) -> Dict[str, float]:
f_plus = np.atleast_1d(job.results.readrkf("Properties", "Fukui Fplus", file="engine")).astype(float)
f_minus = np.atleast_1d(job.results.readrkf("Properties", "Fukui Fminus", file="engine")).astype(float)
f_plus_atom = float(f_plus[atom_index])
f_minus_atom = float(f_minus[atom_index])
return {
"f+": f_plus_atom,
"f-": f_minus_atom,
"f0": 0.5 * (f_plus_atom + f_minus_atom),
"Dual": f_plus_atom - f_minus_atom,
}
rows = []
imgs = {}
for name, job in jobs.items():
row = {"Molecule": name, "Group": substituted_alkenes[name]["group"]}
row.update(read_reactivity_descriptors(job, atom_index=0))
rows.append(row)
imgs[name] = view_atomic_property(
job,
"dual",
direction="along_pca3",
width=400,
padding=1,
label_size=1.8,
colorbar_range=(-0.2, 0.2)
)
df = pd.DataFrame(rows).set_index("Molecule")
df
Group | f+ | f- | f0 | Dual | |
|---|---|---|---|---|---|
Molecule | |||||
Methyl vinyl ether | \(-\mathrm{OCH_3}\) | 0.152974 | 0.238087 | 0.195531 | -0.085113 |
Propene | \(-\mathrm{CH_3}\) | 0.213457 | 0.303950 | 0.258703 | -0.090492 |
Ethene | \(-\mathrm{H}\) | 0.238575 | 0.287503 | 0.263039 | -0.048928 |
Acrylonitrile | \(-\mathrm{CN}\) | 0.267157 | 0.203561 | 0.235359 | 0.063596 |
Acrolein | \(-\mathrm{CHO}\) | 0.251213 | 0.096716 | 0.173965 | 0.154497 |
These can be visualized as part of the wider molecule, for example, the panels below show the dual descriptor in each case:
plot_image_grid(imgs, rows=1);
Compare electron addition and removal at the terminal carbon¶
Both \(f^+\) and \(f^-\) can be appreciable at the same atom. So their difference is therefore especially useful for identifying its dominant character. Negative dual-descriptor values identify a terminal carbon with predominantly nucleophilic character; positive values identify predominantly electrophilic character. In particular, compare the resonance-donating methoxy group with the electron-withdrawing cyano and formyl groups.
def plot_terminal_dual_descriptor(frame: pd.DataFrame) -> Axes:
values = frame["Dual"].to_numpy(dtype=float)
colors = np.where(values < 0.0, "tab:blue", "tab:red")
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(frame["Group"], values, color=colors)
ax.axhline(0.0, color="black", linewidth=0.8)
ax.set_ylabel(r"Dual Descriptor $(f^+-f^-)$")
ax.set_xlabel(r"$\mathrm{R}$ Substituent")
ax.text(0.05, 0.45, "nucleophilic character", color="tab:blue", transform=ax.transAxes, va="top")
ax.text(0.84, 0.35, "electrophilic character", color="tab:red", transform=ax.transAxes, ha="right", va="top")
fig.tight_layout()
return ax
ax = plot_terminal_dual_descriptor(df)
ax;
Interpretation and limitations¶
For methyl vinyl ether, resonance donation from oxygen gives the terminal carbon nucleophilic character: it donates electron density to an attacking electrophile. The methyl group in propene also gives the terminal carbon nucleophilic character through electron donation by hyperconjugation. Ethene provides the unsubstituted reference and has a weaker negative dual descriptor at its carbon atoms. By contrast, the electron-withdrawing cyano group in acrylonitrile reverses the response and makes the terminal carbon electrophilic, so it accepts electron density from an attacking nucleophile. The formyl group in acrolein produces the strongest electrophilic response of this series at the terminal carbon, consistent with its role as the beta carbon in conjugate nucleophilic addition.
Condensed Fukui functions are normalized within each molecule and are best used to compare sites or qualitative response patterns. They do not directly predict reaction rates. Their numerical values also depend on the atomic-charge model; here they reflect AIMNet2’s learned charges rather than a particular quantum-chemical population analysis.
See also¶
References¶
Python Script¶
#!/usr/bin/env python
# coding: utf-8
# ## Predicting Reactivity via Condensed Fukui Functions
#
# This example uses the AIMNet2-NSE machine learning potential to compare a series of molecules with the common vinyl fragment $\mathrm{CH_2{=}CH{-}R}$. The substituent $R$ ranges from electron donating to electron withdrawing. This will demonstrate how the terminal vinyl carbon changes from being predominantly nucleophilic to predominantly electrophilic in nature.
#
# For an atom $A$, the condensed Fukui functions describe its response to adding or removing an electron. $f_A^+$ measures the response to electron addition and is large at electrophilic sites, which accept electron density from nucleophiles. $f_A^-$ measures the response to electron removal and is large at nucleophilic sites, which donate electron density to electrophiles.
#
# The dual descriptor is $\Delta f_A=f_A^+-f_A^-$. Positive values indicate predominantly electrophilic character, while negative values indicate predominantly nucleophilic character.
#
# AIMNet2 evaluates the electron-added and electron-removed states at the optimized neutral geometry. These same calculations provide the vertical electron affinity and ionization energy.
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.axes import Axes
from scm.base import ChemicalSystem, Units
from scm.plams import AMSJob, Settings, plot_image_grid, view, view_atomic_property
# ### Build the substituted alkenes
#
# The SMILES strings are written from the terminal vinyl carbon, making it atom 1 in every molecule. Methoxy is a resonance donor, methyl is a weak donor through hyperconjugation, and cyano and formyl are electron-withdrawing groups. Ethene provides the unsubstituted reference.
substituted_alkenes = {
"Methyl vinyl ether": {"smiles": "C=COC", "group": r"$-\mathrm{OCH_3}$"},
"Propene": {"smiles": "C=CC", "group": r"$-\mathrm{CH_3}$"},
"Ethene": {"smiles": "C=C", "group": r"$-\mathrm{H}$"},
"Acrylonitrile": {"smiles": "C=CC#N", "group": r"$-\mathrm{CN}$"},
"Acrolein": {"smiles": "C=CC=O", "group": r"$-\mathrm{CHO}$"},
}
molecules = {
name: ChemicalSystem.from_smiles(details["smiles"])
for name, details in substituted_alkenes.items()
}
plot_image_grid(
{name: view(molecule, direction="along_pca3") for name, molecule in molecules.items()}, rows=1
, save_path="picture1.png");
# ### Optimize the neutral structures and calculate the descriptors
#
# When AMS Properties%Other is requested, the condensed Fukui functions are automatically calculated for an aperiodic ML potential that accepts total charge and returns individual atomic charges. This automatic behavior can be controlled with `Fukui%Enabled`. During a geometry optimization, the charged states are evaluated only after optimizing the neutral system. Their geometry is held fixed, so all reported energy differences are vertical.
def reactivity_settings() -> Settings:
settings = Settings()
settings.input.ams.Task = "GeometryOptimization"
settings.input.MLPotential.Model = "AIMNet2-NSE"
settings.runscript.nproc = 1
return settings
jobs = {}
for name, molecule in molecules.items():
job_name = name.lower().replace(" ", "_")
jobs[name] = AMSJob(name=job_name, molecule=molecule, settings=reactivity_settings())
jobs[name].run();
# ### Read and plot descriptor results
#
# The Fukui arrays are stored in the MLPotential engine results. The helper below accepts a zero-based Python atom index. Thus, index 0 corresponds to atom 1: the terminal $\mathrm{CH_2}$ carbon in every structure.
def read_reactivity_descriptors(job: AMSJob, atom_index: int) -> Dict[str, float]:
f_plus = np.atleast_1d(job.results.readrkf("Properties", "Fukui Fplus", file="engine")).astype(float)
f_minus = np.atleast_1d(job.results.readrkf("Properties", "Fukui Fminus", file="engine")).astype(float)
f_plus_atom = float(f_plus[atom_index])
f_minus_atom = float(f_minus[atom_index])
return {
"f+": f_plus_atom,
"f-": f_minus_atom,
"f0": 0.5 * (f_plus_atom + f_minus_atom),
"Dual": f_plus_atom - f_minus_atom,
}
rows = []
imgs = {}
for name, job in jobs.items():
row = {"Molecule": name, "Group": substituted_alkenes[name]["group"]}
row.update(read_reactivity_descriptors(job, atom_index=0))
rows.append(row)
imgs[name] = view_atomic_property(
job,
"dual",
direction="along_pca3",
width=400,
padding=1,
label_size=1.8,
colorbar_range=(-0.2, 0.2)
)
df = pd.DataFrame(rows).set_index("Molecule")
print(df)
# These can be visualized as part of the wider molecule, for example, the panels below show the dual descriptor in each case:
plot_image_grid(imgs, rows=1, save_path="picture2.png");
# ### Compare electron addition and removal at the terminal carbon
#
# Both $f^+$ and $f^-$ can be appreciable at the same atom. So their difference is therefore especially useful for identifying its dominant character. Negative dual-descriptor values identify a terminal carbon with predominantly nucleophilic character; positive values identify predominantly electrophilic character. In particular, compare the resonance-donating methoxy group with the electron-withdrawing cyano and formyl groups.
def plot_terminal_dual_descriptor(frame: pd.DataFrame) -> Axes:
values = frame["Dual"].to_numpy(dtype=float)
colors = np.where(values < 0.0, "tab:blue", "tab:red")
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(frame["Group"], values, color=colors)
ax.axhline(0.0, color="black", linewidth=0.8)
ax.set_ylabel(r"Dual Descriptor $(f^+-f^-)$")
ax.set_xlabel(r"$\mathrm{R}$ Substituent")
ax.text(0.05, 0.45, "nucleophilic character", color="tab:blue", transform=ax.transAxes, va="top")
ax.text(0.84, 0.35, "electrophilic character", color="tab:red", transform=ax.transAxes, ha="right", va="top")
fig.tight_layout()
return ax
ax = plot_terminal_dual_descriptor(df)
ax;
ax.figure.savefig("picture3.png")
# ### Interpretation and limitations
#
# For methyl vinyl ether, resonance donation from oxygen gives the terminal carbon nucleophilic character: it donates electron density to an attacking electrophile. The methyl group in propene also gives the terminal carbon nucleophilic character through electron donation by hyperconjugation. Ethene provides the unsubstituted reference and has a weaker negative dual descriptor at its carbon atoms. By contrast, the electron-withdrawing cyano group in acrylonitrile reverses the response and makes the terminal carbon electrophilic, so it accepts electron density from an attacking nucleophile. The formyl group in acrolein produces the strongest electrophilic response of this series at the terminal carbon, consistent with its role as the beta carbon in conjugate nucleophilic addition.
#
# Condensed Fukui functions are normalized within each molecule and are best used to compare sites or qualitative response patterns. They do not directly predict reaction rates. Their numerical values also depend on the atomic-charge model; here they reflect AIMNet2's learned charges rather than a particular quantum-chemical population analysis.