AMS transition state workflow

Example illustrating a Diels-Alder addition between cyclopentadiene and acrylonitrile.

In this example, a preliminary biased UFF MD simulation is run to orient the two reactant molecules near the transition state.

Then a transition state search is performed.

The PESExploration LandscapeRefinement is used to find the two minima on either side of the transition state.

See also

The AMS biased MD / PLUMED example is similar, but in that example the reaction happens during the MD simulation. That type of simulation can be performed when the potential is reactive.

In the current example, the MD simulation is performed with UFF which is not reactive.

See also

Example: Molecule substitution: Attach ligands to substrates

To follow along, either

• Download diels_alder_addition.py (run as $AMSBIN/amspython diels_alder_addition.py).

• Download diels_alder_addition.ipynb (see also: how to install Jupyterlab in AMS)

Initial imports

from scm.plams import *
import numpy as np
import os
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 uses the PESExploration LandscapeRefinement tool to get the two corresponding minima and energy landscape. Alternatively, one could also do an IRC (intrinsic reaction coordinate) calculation.

def addition(
    mol1: Molecule,
    mol2: Molecule,
    covalent_radius_multiplier: float = 1.2,
):
    mol1_ind = get_active_i(mol1)
    mol2_ind = get_active_i(mol2)
    if len(mol1_ind) != 2:
        raise ValueError(
            f"Set at.properties.active_bond for two atoms in mol1. Current mol1_ind: {mol1_ind}"
        )
    if len(mol2_ind) != 2:
        raise ValueError(
            f"Set at.properties.active_bond for two atoms in mol2. Current mol2_ind: {mol2_ind}"
        )

    target_distances = []

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

    combined = mol1.copy()
    combined.add_molecule(mol2, margin=2)
    combined.properties.charge = mol1.properties.get("charge", 0) + mol2.properties.get("charge", 0)

    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,
    )
    mol = preliminary_md_results.get_main_molecule()

    ts_search_results = ts_search(mol, mol1_ind, mol2_ind)
    mol = ts_search_results.get_main_molecule()

    relax_from_saddle_results = relax_from_saddle(mol)

    return preliminary_md_results, ts_search_results, relax_from_saddle_results


def ts_search(
    molecule, atom_indices_1: List[int], atom_indices_2: List[int], settings:Settings=None
) -> AMSResults:
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        # settings.runscript.nproc = 1
        settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"
        settings.input.ams.Properties.NormalModes = "Yes"
        settings.input.ams.GeometryOptimization.Convergence.Quality = "Good"

    s = settings.copy()
    s.input.ams.task = "TransitionStateSearch"
    s.input.ams.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 relax_from_saddle(molecule:Molecule, settings:Settings=None) -> AMSResults:
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"
        # settings.runscript.nproc = 1

    s = settings.copy()
    s.input.ams.task = "PESExploration"
    s.input.ams.PESExploration.Job = "LandscapeRefinement"
    s.input.ams.PESExploration.LandscapeRefinement.RelaxFromSaddlePoint = "T"

    dummy = from_smiles("[HH]")
    m = {"state1 ts=Yes": molecule, "": dummy}

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

    return job.results


def irc(molecule:Molecule, settings:Settings=None):
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        # settings.runscript.nproc = 1

    s = settings.copy()
    s.input.ams.task = "IRC"
    s.input.ams.IRC.MinEnergyProfile = "Yes"

    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: Molecule,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    target_distances: List[float],
    nsteps: int = 10000,
    kappa: float = 100000,
    settings: Settings = None,
    temperature: float = 300,
) -> Molecule:
    if settings is None:
        settings = Settings()
        settings.input.ForceField.Type = "UFF"
        settings.runscript.nproc = 1

    plumed_input = "\n"
    for a1, a2, d in zip(atom_indices_1, atom_indices_2, target_distances):
        # current_d = molecule[a1].distance_to(molecule[a2])
        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"

    plumed_input += "   End"
    settings.input.ams.MolecularDynamics.Plumed.Input = plumed_input

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

    return job.results


def set_active(mol: Molecule, indices: List[int]):
    for at in mol:
        if "active_bond" in at.properties:
            del at.properties["active_bond"]
    for i, ind in enumerate(indices, 1):
        mol[ind].properties.active_bond = i


def get_active_i(mol: Molecule) -> List[int]:
    d = {}
    for i, at in enumerate(mol, 1):
        if 'active_bond' in at.properties and at.properties.active_bond:
            d[i] = at.properties.active_bond

    return sorted(d, key=lambda x: d[x])


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

Run the calculations

init()

PLAMS working folder: AMSTSWorkflow/plams_workdir

diene_smiles = "C1C=CC=C1"
diene = from_smiles(diene_smiles)  # carbon 2, 4 will form bonds. Do diene.write('diene.xyz') and open diene.xyz in the AMS GUI to find out which atom indices are correct.
set_active(diene, [2, 4])

plot_molecule(diene)
plt.title('Diene (cyclopentadiene)');

dienophile = from_smiles("N#CC=C")  # carbon 1, 2 will form bonds
set_active(dienophile, [1, 2])

plot_molecule(dienophile)
plt.title('Dienophile (acrylonitrile)');

preliminary_md_results, ts_search_results, relax_from_saddle_results = addition(diene, dienophile)

[28.03|14:41:11] JOB preliminary_md STARTED
[28.03|14:41:11] JOB preliminary_md RUNNING
[28.03|14:41:18] JOB preliminary_md FINISHED
[28.03|14:41:18] JOB preliminary_md SUCCESSFUL
[28.03|14:41:18] JOB ts_search STARTED
[28.03|14:41:18] JOB ts_search RUNNING
[28.03|14:41:21] JOB ts_search FINISHED
[28.03|14:41:21] JOB ts_search SUCCESSFUL
[28.03|14:41:21] JOB refinement STARTED
[28.03|14:41:21] JOB refinement RUNNING
[28.03|14:41:27] JOB refinement FINISHED
[28.03|14:41:27] JOB refinement SUCCESSFUL

Preliminary biased MD results (UFF)

final_md_system = preliminary_md_results.get_main_molecule()
plot_molecule(final_md_system)
plt.title('Final system from preliminary biased MD');

TS search results (DFTB)

final_ts_system = ts_search_results.get_main_molecule()
plot_molecule(final_ts_system)
plt.title("DFTB-optimized transition state");

Energy landscape refinement results (DFTB)

landscape = relax_from_saddle_results.get_energy_landscape()
print(landscape)

All stationary points:
======================
State 1: C8H9N local minimum @ -24.90376082 Hartree (found 1 times, results on State-1_MIN)
State 2: C8H9N local minimum @ -24.83422501 Hartree (found 1 times, results on State-2_MIN)
State 3: C8H9N transition state @ -24.82034339 Hartree (found 1 times, results on State-3_TS_1-2)
  +- Reactants: State 1: C8H9N local minimum @ -24.90376082 Hartree (found 1 times, results on State-1_MIN)
     Products:  State 2: C8H9N local minimum @ -24.83422501 Hartree (found 1 times, results on State-2_MIN)
     Prefactors: 0.000E+00:0.000E+00 s^-1
     Barriers: 2.270:0.378 eV

Above we see that the forward and backward barriers are 2.27 and 0.39 eV, respectively.

Ha2eV = Units.convert(1.0, 'hartree', 'eV')
energies = landscape[1].energy, landscape[3].energy, landscape[2].energy
energies = (np.array(energies) - landscape[1].energy) * Ha2eV
plt.plot(energies)
plt.ylabel('Relative energy (eV)')
plt.xticks([0, 1, 2], ['State 1 (min)', 'State 3 (TS)', 'State 2 (min)']);

plot_molecule(landscape[1].molecule)
plt.title('State 1 (minimum)');

plot_molecule(landscape[2].molecule)
plt.title('State 2 (minimum)');

plot_molecule(landscape[3].molecule)
plt.title('State 3 (Transition state)');

Finish plams

finish()

[28.03|14:41:28] PLAMS run finished. Goodbye

Complete Python code

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

# ## Initial imports

from typing import List

import matplotlib.pyplot as plt
import numpy as np
from scm.plams import *

# ## 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 uses the PESExploration LandscapeRefinement tool to get the two corresponding minima and energy landscape. Alternatively, one could also do an IRC (intrinsic reaction coordinate) calculation.


def addition(
    mol1: Molecule,
    mol2: Molecule,
    covalent_radius_multiplier: float = 1.2,
):
    mol1_ind = get_active_i(mol1)
    mol2_ind = get_active_i(mol2)
    if len(mol1_ind) != 2:
        raise ValueError(f"Set at.properties.active_bond for two atoms in mol1. Current mol1_ind: {mol1_ind}")
    if len(mol2_ind) != 2:
        raise ValueError(f"Set at.properties.active_bond for two atoms in mol2. Current mol2_ind: {mol2_ind}")

    target_distances = []

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

    combined = mol1.copy()
    combined.add_molecule(mol2, margin=2)
    combined.properties.charge = mol1.properties.get("charge", 0) + mol2.properties.get("charge", 0)

    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,
    )
    mol = preliminary_md_results.get_main_molecule()

    ts_search_results = ts_search(mol, mol1_ind, mol2_ind)
    mol = ts_search_results.get_main_molecule()

    relax_from_saddle_results = relax_from_saddle(mol)

    return preliminary_md_results, ts_search_results, relax_from_saddle_results


def ts_search(molecule, atom_indices_1: List[int], atom_indices_2: List[int], settings: Settings = None) -> AMSResults:
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        # settings.runscript.nproc = 1
        settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"
        settings.input.ams.Properties.NormalModes = "Yes"
        settings.input.ams.GeometryOptimization.Convergence.Quality = "Good"

    s = settings.copy()
    s.input.ams.task = "TransitionStateSearch"
    s.input.ams.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 relax_from_saddle(molecule: Molecule, settings: Settings = None) -> AMSResults:
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"
        # settings.runscript.nproc = 1

    s = settings.copy()
    s.input.ams.task = "PESExploration"
    s.input.ams.PESExploration.Job = "LandscapeRefinement"
    s.input.ams.PESExploration.LandscapeRefinement.RelaxFromSaddlePoint = "T"

    dummy = from_smiles("[HH]")
    m = {"state1 ts=Yes": molecule, "": dummy}

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

    return job.results


def irc(molecule: Molecule, settings: Settings = None):
    if settings is None:
        settings = Settings()
        settings.input.DFTB
        # settings.runscript.nproc = 1

    s = settings.copy()
    s.input.ams.task = "IRC"
    s.input.ams.IRC.MinEnergyProfile = "Yes"

    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: Molecule,
    atom_indices_1: List[int],
    atom_indices_2: List[int],
    target_distances: List[float],
    nsteps: int = 10000,
    kappa: float = 100000,
    settings: Settings = None,
    temperature: float = 300,
) -> Molecule:
    if settings is None:
        settings = Settings()
        settings.input.ForceField.Type = "UFF"
        settings.runscript.nproc = 1

    plumed_input = "\n"
    for a1, a2, d in zip(atom_indices_1, atom_indices_2, target_distances):
        # current_d = molecule[a1].distance_to(molecule[a2])
        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 += "\n"

    plumed_input += "   End"
    settings.input.ams.MolecularDynamics.Plumed.Input = plumed_input

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

    return job.results


def set_active(mol: Molecule, indices: List[int]):
    for at in mol:
        if "active_bond" in at.properties:
            del at.properties["active_bond"]
    for i, ind in enumerate(indices, 1):
        mol[ind].properties.active_bond = i


def get_active_i(mol: Molecule) -> List[int]:
    d = {}
    for i, at in enumerate(mol, 1):
        if "active_bond" in at.properties and at.properties.active_bond:
            d[i] = at.properties.active_bond

    return sorted(d, key=lambda x: d[x])


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


# ## Run the calculations

init()


diene_smiles = "C1C=CC=C1"
diene = from_smiles(
    diene_smiles
)  # carbon 2, 4 will form bonds. Do diene.write('diene.xyz') and open diene.xyz in the AMS GUI to find out which atom indices are correct.
set_active(diene, [2, 4])


plot_molecule(diene)
plt.title("Diene (cyclopentadiene)")


dienophile = from_smiles("N#CC=C")  # carbon 1, 2 will form bonds
set_active(dienophile, [1, 2])


plot_molecule(dienophile)
plt.title("Dienophile (acrylonitrile)")


preliminary_md_results, ts_search_results, relax_from_saddle_results = addition(diene, dienophile)


# ## Preliminary biased MD results (UFF)

final_md_system = preliminary_md_results.get_main_molecule()
plot_molecule(final_md_system)
plt.title("Final system from preliminary biased MD")


# ## TS search results (DFTB)

final_ts_system = ts_search_results.get_main_molecule()
plot_molecule(final_ts_system)
plt.title("DFTB-optimized transition state")


# ## Energy landscape refinement results (DFTB)

landscape = relax_from_saddle_results.get_energy_landscape()
print(landscape)


# Above we see that the forward and backward barriers are 2.27 and 0.39 eV, respectively.

Ha2eV = Units.convert(1.0, "hartree", "eV")
energies = landscape[1].energy, landscape[3].energy, landscape[2].energy
energies = (np.array(energies) - landscape[1].energy) * Ha2eV
plt.plot(energies)
plt.ylabel("Relative energy (eV)")
plt.xticks([0, 1, 2], ["State 1 (min)", "State 3 (TS)", "State 2 (min)"])


plot_molecule(landscape[1].molecule)
plt.title("State 1 (minimum)")


plot_molecule(landscape[2].molecule)
plt.title("State 2 (minimum)")


plot_molecule(landscape[3].molecule)
plt.title("State 3 (Transition state)")


# ## Finish plams

finish()