Automated Transition-State Workflow for a Diels-Alder Addition¶

This example illustrates a Diels-Alder addition between cyclopentadiene and acrylonitrile.

A preliminary biased UFF molecular-dynamics simulation is used to orient the two reactant molecules near the transition state, after which a transition-state search is performed. The PESExploration LandscapeRefinement step is then used to locate the two minima on either side of the transition state.

A related workflow is shown in PLUMED Biasing in AMS Molecular Dynamics, where the reaction is driven directly during the MD simulation with a reactive potential.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

../_images/ams_ts_workflow_15_0_67e5b9ba.png

Initial imports¶

from scm.base import Atom, ChemicalSystem
from scm.inputs import AMS, DFTB, ForceField
from scm.plams import AMSJob, AMSNVTJob, AMSResults, PeriodicTable, Settings, Units, plot_image_grid, view
import numpy as np
import matplotlib.pyplot as plt
from typing import List

Function definitions¶

The addition() function

performs a preliminary MD simulation with UFF to arrange the two reactant molecules at an approximate transition state,
then does a transition state search with DFTB, specifying the known reaction coordinate (bond formation),
then follows the transition state downhill in both directions with an intrinsic reaction coordinate (IRC) calculation to obtain the adjacent minima and energy profile.

def addition(
    mol1: ChemicalSystem,
    mol2: ChemicalSystem,
    mol1_ind: List[int],
    mol2_ind: List[int],
    covalent_radius_multiplier: float = 1.2,
):
    if len(mol1_ind) != 2:
        raise ValueError(
            f"Provide two forming-bond atom indices for mol1. Current mol1_ind: {mol1_ind}"
        )
    if len(mol2_ind) != 2:
        raise ValueError(
            f"Provide two forming-bond atom indices for mol2. Current mol2_ind: {mol2_ind}"
        )

    target_distances = []

    for i1, i2 in zip(mol1_ind, mol2_ind):
        r = get_atom_radius(mol1.atoms[i1 - 1]) + get_atom_radius(mol2.atoms[i2 - 1])
        target_distances.append(r * covalent_radius_multiplier)

    combined = mol1.copy()
    translated_mol2 = mol2.copy()
    x_offset = max(atom.coords[0] for atom in combined.atoms) - min(
        atom.coords[0] for atom in translated_mol2.atoms
    ) + 2.0
    translated_mol2.translate((x_offset, 0.0, 0.0))
    combined += translated_mol2

    mol2_ind = [x + len(mol1) for x in mol2_ind]

    preliminary_md_results = preliminary_md(
        combined,
        mol1_ind,
        mol2_ind,
        target_distances=target_distances,
        temperature=600,
        nsteps=20000,
    )
    system = preliminary_md_results.get_main_system()

    ts_search_results = ts_search(system, mol1_ind, mol2_ind)
    system = ts_search_results.get_main_system()

    irc_results = irc(system)

    return preliminary_md_results, ts_search_results, irc_results


def ts_search(
    molecule: ChemicalSystem,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    settings: AMS = None,
) -> AMSResults:
    if settings is None:
        settings = AMS(Engine=DFTB())
        settings.GeometryOptimization.InitialHessian.Type = "Calculate"
        settings.Properties.NormalModes = True
        settings.GeometryOptimization.Convergence.Quality = "Good"

    s = settings.model_copy(deep=True)
    s.Task = "TransitionStateSearch"
    s.TransitionStateSearch.ReactionCoordinate.Distance = [
        f"{a1} {a2} 1.0" for a1, a2 in zip(atom_indices_1, atom_indices_2)
    ]

    job = AMSJob(settings=s, name="ts_search", molecule=molecule)
    job.run()

    return job.results


def irc(molecule: ChemicalSystem, settings: AMS = None):
    if settings is None:
        settings = AMS(Engine=DFTB())

    s = settings.model_copy(deep=True)
    s.Task = "IRC"
    s.IRC.MinEnergyProfile = True

    job = AMSJob(settings=s, name="irc", molecule=molecule)
    job.run()

    converged = job.results.get_history_property("Converged")
    direction = job.results.get_history_property("IRCDirection")
    energies = job.results.get_history_property("Energy")
    ind = dict()
    energy = dict()
    for i, (c, d, e) in enumerate(zip(converged, direction, energies)):
        if c:
            ind[d] = i + 1
            energy[d] = e

    min1 = job.results.get_history_molecule(ind[1])
    min2 = job.results.get_history_molecule(ind[2])
    ts = job.results.get_input_molecule()

    return min1, energy[1], min2, energy[2], ts, energies[0]


def preliminary_md(
    molecule: ChemicalSystem,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    target_distances: List[float],
    nsteps: int = 10000,
    kappa: float = 100000,
    settings: AMS = None,
    temperature: float = 300,
) -> AMSResults:
    use_serial = settings is None
    if settings is None:
        settings = AMS(Engine=ForceField(Type="UFF"))
    settings = settings.model_copy(deep=True)

    plumed_input = "\n"
    for a1, a2, d in zip(atom_indices_1, atom_indices_2, target_distances):
        plumed_input += f"DISTANCE ATOMS={a1},{a2} LABEL=d_{a1}_{a2}\n"
        plumed_input += f"MOVINGRESTRAINT ARG=d_{a1}_{a2}"
        plumed_input += f" STEP0=1 AT0={d * 0.1} KAPPA0=0"
        plumed_input += f" STEP1={1 * nsteps // 4} KAPPA1={kappa / 1000}"
        plumed_input += f" STEP2={2 * nsteps // 4} KAPPA2={kappa / 100}"
        plumed_input += f" STEP3={3 * nsteps // 4} KAPPA3={kappa / 10}"
        plumed_input += f" STEP4={4 * nsteps // 4} KAPPA4={kappa}"
        plumed_input += f"\n"

    settings.MolecularDynamics.Plumed.Input = plumed_input
    job_settings = Settings()
    job_settings.input = settings
    if use_serial:
        job_settings.runscript.nproc = 1

    job = AMSNVTJob(
        name="preliminary_md",
        settings=job_settings,
        molecule=molecule,
        nsteps=nsteps,
        temperature=3 * [temperature] + [1],
    )
    job.run()

    return job.results


def get_atom_radius(at: Atom) -> float:
    return PeriodicTable.get_radius(at.symbol)

Run the calculations¶

diene_smiles = "C1C=CC=C1"
diene = ChemicalSystem.from_smiles(diene_smiles)
diene_indices = [2, 5]  # terminal carbon atoms of the conjugated diene

dienophile = ChemicalSystem.from_smiles("N#CC=C")
dienophile_indices = [4, 3]  # alkene carbon atoms, terminal carbon first

plot_image_grid(
    {"Diene (cyclopentadiene)": view(diene), "Dienophile (acrylonitrile)": view(dienophile)}
);

preliminary_md_results, ts_search_results, irc_results = addition(
    diene, dienophile, diene_indices, dienophile_indices
)

[16.10|09:20:22] JOB preliminary_md STARTED
[16.10|09:20:22] JOB preliminary_md RUNNING
[16.10|09:20:28] JOB preliminary_md FINISHED
[16.10|09:20:28] JOB preliminary_md SUCCESSFUL
[16.10|09:20:28] JOB ts_search STARTED
[16.10|09:20:28] JOB ts_search RUNNING
[16.10|09:20:40] JOB ts_search FINISHED
[16.10|09:20:40] JOB ts_search SUCCESSFUL
[16.10|09:20:40] JOB irc STARTED
[16.10|09:20:40] JOB irc RUNNING
[16.10|09:20:54] JOB irc FINISHED
[16.10|09:20:54] JOB irc SUCCESSFUL

Preliminary biased MD (UFF) & TS search (DFTB) results¶

final_md_system = preliminary_md_results.get_main_system()
final_ts_system = ts_search_results.get_main_system()
plot_image_grid(
    {
        "Final system from preliminary biased MD": view(final_md_system),
        "DFTB-optimized transition state": view(final_ts_system),
    }
);

Intrinsic reaction coordinate results (DFTB)¶

minimum_1, energy_1, minimum_2, energy_2, transition_state, transition_state_energy = irc_results

The IRC follows the transition state downhill in both directions to the adjacent minima.

hartree_to_ev = Units.convert(1.0, "hartree", "eV")
energies = np.array([energy_1, transition_state_energy, energy_2])
energies = (energies - min(energy_1, energy_2)) * hartree_to_ev
fig, ax = plt.subplots()
ax.plot(energies)
ax.set_ylabel("Relative energy (eV)")
ax.set_xticks([0, 1, 2], ["Minimum 1", "Transition state", "Minimum 2"])
ax;

plot_image_grid(
    {
        "Minimum 1": view(minimum_1, guess_bonds=True),
        "Transition state": view(transition_state, guess_bonds=True),
        "Minimum 2": view(minimum_2, guess_bonds=True),
    },
    rows=1,
);

Python Script¶

#!/usr/bin/env python
# coding: utf-8

# ## Initial imports

from scm.base import Atom, ChemicalSystem
from scm.inputs import AMS, DFTB, ForceField
from scm.plams import AMSJob, AMSNVTJob, AMSResults, PeriodicTable, Settings, Units, plot_image_grid, view
import numpy as np
import matplotlib.pyplot as plt
from typing import List


# ## Function definitions
# 
# The ``addition()`` function
# 
# * performs a preliminary MD simulation with UFF to arrange the two reactant molecules at an approximate transition state,
# 
# * then does a transition state search with DFTB, specifying the known reaction coordinate (bond formation),
# 
# * then follows the transition state downhill in both directions with an intrinsic reaction coordinate (IRC) calculation to obtain the adjacent minima and energy profile.

def addition(
    mol1: ChemicalSystem,
    mol2: ChemicalSystem,
    mol1_ind: List[int],
    mol2_ind: List[int],
    covalent_radius_multiplier: float = 1.2,
):
    if len(mol1_ind) != 2:
        raise ValueError(
            f"Provide two forming-bond atom indices for mol1. Current mol1_ind: {mol1_ind}"
        )
    if len(mol2_ind) != 2:
        raise ValueError(
            f"Provide two forming-bond atom indices for mol2. Current mol2_ind: {mol2_ind}"
        )

    target_distances = []

    for i1, i2 in zip(mol1_ind, mol2_ind):
        r = get_atom_radius(mol1.atoms[i1 - 1]) + get_atom_radius(mol2.atoms[i2 - 1])
        target_distances.append(r * covalent_radius_multiplier)

    combined = mol1.copy()
    translated_mol2 = mol2.copy()
    x_offset = max(atom.coords[0] for atom in combined.atoms) - min(
        atom.coords[0] for atom in translated_mol2.atoms
    ) + 2.0
    translated_mol2.translate((x_offset, 0.0, 0.0))
    combined += translated_mol2

    mol2_ind = [x + len(mol1) for x in mol2_ind]

    preliminary_md_results = preliminary_md(
        combined,
        mol1_ind,
        mol2_ind,
        target_distances=target_distances,
        temperature=600,
        nsteps=20000,
    )
    system = preliminary_md_results.get_main_system()

    ts_search_results = ts_search(system, mol1_ind, mol2_ind)
    system = ts_search_results.get_main_system()

    irc_results = irc(system)

    return preliminary_md_results, ts_search_results, irc_results


def ts_search(
    molecule: ChemicalSystem,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    settings: AMS = None,
) -> AMSResults:
    if settings is None:
        settings = AMS(Engine=DFTB())
        settings.GeometryOptimization.InitialHessian.Type = "Calculate"
        settings.Properties.NormalModes = True
        settings.GeometryOptimization.Convergence.Quality = "Good"

    s = settings.model_copy(deep=True)
    s.Task = "TransitionStateSearch"
    s.TransitionStateSearch.ReactionCoordinate.Distance = [
        f"{a1} {a2} 1.0" for a1, a2 in zip(atom_indices_1, atom_indices_2)
    ]

    job = AMSJob(settings=s, name="ts_search", molecule=molecule)
    job.run()

    return job.results


def irc(molecule: ChemicalSystem, settings: AMS = None):
    if settings is None:
        settings = AMS(Engine=DFTB())

    s = settings.model_copy(deep=True)
    s.Task = "IRC"
    s.IRC.MinEnergyProfile = True

    job = AMSJob(settings=s, name="irc", molecule=molecule)
    job.run()

    converged = job.results.get_history_property("Converged")
    direction = job.results.get_history_property("IRCDirection")
    energies = job.results.get_history_property("Energy")
    ind = dict()
    energy = dict()
    for i, (c, d, e) in enumerate(zip(converged, direction, energies)):
        if c:
            ind[d] = i + 1
            energy[d] = e

    min1 = job.results.get_history_molecule(ind[1])
    min2 = job.results.get_history_molecule(ind[2])
    ts = job.results.get_input_molecule()

    return min1, energy[1], min2, energy[2], ts, energies[0]


def preliminary_md(
    molecule: ChemicalSystem,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    target_distances: List[float],
    nsteps: int = 10000,
    kappa: float = 100000,
    settings: AMS = None,
    temperature: float = 300,
) -> AMSResults:
    use_serial = settings is None
    if settings is None:
        settings = AMS(Engine=ForceField(Type="UFF"))
    settings = settings.model_copy(deep=True)

    plumed_input = "\n"
    for a1, a2, d in zip(atom_indices_1, atom_indices_2, target_distances):
        plumed_input += f"DISTANCE ATOMS={a1},{a2} LABEL=d_{a1}_{a2}\n"
        plumed_input += f"MOVINGRESTRAINT ARG=d_{a1}_{a2}"
        plumed_input += f" STEP0=1 AT0={d * 0.1} KAPPA0=0"
        plumed_input += f" STEP1={1 * nsteps // 4} KAPPA1={kappa / 1000}"
        plumed_input += f" STEP2={2 * nsteps // 4} KAPPA2={kappa / 100}"
        plumed_input += f" STEP3={3 * nsteps // 4} KAPPA3={kappa / 10}"
        plumed_input += f" STEP4={4 * nsteps // 4} KAPPA4={kappa}"
        plumed_input += f"\n"

    settings.MolecularDynamics.Plumed.Input = plumed_input
    job_settings = Settings()
    job_settings.input = settings
    if use_serial:
        job_settings.runscript.nproc = 1

    job = AMSNVTJob(
        name="preliminary_md",
        settings=job_settings,
        molecule=molecule,
        nsteps=nsteps,
        temperature=3 * [temperature] + [1],
    )
    job.run()

    return job.results


def get_atom_radius(at: Atom) -> float:
    return PeriodicTable.get_radius(at.symbol)


# ## Run the calculations

diene_smiles = "C1C=CC=C1"
diene = ChemicalSystem.from_smiles(diene_smiles)
diene_indices = [2, 5]  # terminal carbon atoms of the conjugated diene


dienophile = ChemicalSystem.from_smiles("N#CC=C")
dienophile_indices = [4, 3]  # alkene carbon atoms, terminal carbon first


plot_image_grid(
    {"Diene (cyclopentadiene)": view(diene), "Dienophile (acrylonitrile)": view(dienophile)}
, save_path="picture1.png");


preliminary_md_results, ts_search_results, irc_results = addition(
    diene, dienophile, diene_indices, dienophile_indices
)


# ## Preliminary biased MD (UFF) & TS search (DFTB) results

final_md_system = preliminary_md_results.get_main_system()
final_ts_system = ts_search_results.get_main_system()
plot_image_grid(
    {
        "Final system from preliminary biased MD": view(final_md_system),
        "DFTB-optimized transition state": view(final_ts_system),
    }
, save_path="picture2.png");


# ## Intrinsic reaction coordinate results (DFTB)

minimum_1, energy_1, minimum_2, energy_2, transition_state, transition_state_energy = irc_results


# The IRC follows the transition state downhill in both directions to the adjacent minima.

hartree_to_ev = Units.convert(1.0, "hartree", "eV")
energies = np.array([energy_1, transition_state_energy, energy_2])
energies = (energies - min(energy_1, energy_2)) * hartree_to_ev
fig, ax = plt.subplots()
ax.plot(energies)
ax.set_ylabel("Relative energy (eV)")
ax.set_xticks([0, 1, 2], ["Minimum 1", "Transition state", "Minimum 2"])
ax;
ax.figure.savefig("picture3.png")


plot_image_grid(
    {
        "Minimum 1": view(minimum_1, guess_bonds=True),
        "Transition state": view(transition_state, guess_bonds=True),
        "Minimum 2": view(minimum_2, guess_bonds=True),
    },
    rows=1,
 save_path="picture4.png");

Automated Transition-State Workflow for a Diels-Alder Addition¶

Initial imports¶

Function definitions¶

Run the calculations¶

Preliminary biased MD (UFF) & TS search (DFTB) results¶

Intrinsic reaction coordinate results (DFTB)¶

See also¶

Related examples¶

Related tutorials¶

Related documentation¶

Python Script¶