Transition state search for polymerization: Ziegler Natta catalyst¶

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

# ## Plan
#
# This tutorial shows how to
#
# * optimize the reactant state,
# * perform a PES Scan to get a good starting structure for transition state search,
# * perform a transition state search,
# * extract the barrier
#
# Note that this tutorial only shows one half of a reaction.

# ## Initial imports

from scm.libbase import ChemicalSystem
from scm.plams import view, Settings, AMSJob


# ## Initial system

mol = ChemicalSystem(
    """
System
    Charge 1.0
    Atoms
        Zr       1.60682759      -9.09697636       1.29667400
        C       2.89496678     -10.06750666       3.11673910
        C       2.84805610      -9.90469975      -0.63831357
        C       1.91952398     -11.03393562       2.76898899
        C       2.21720521      -8.93716871       3.64709384
        C       2.05919385      -8.80831375      -1.07706558
        C       1.96898612     -10.95146163      -0.26897070
        C       0.63990901     -10.49043590       3.04466191
        C       0.82713373      -9.20286661       3.61113029
        C       0.69607878      -9.17611676      -0.96529901
        C       0.63812917     -10.49403820      -0.44270184
        H       3.96827494     -10.23973626       3.15013890
        H       3.92582990      -9.99794535      -0.74784951
        H       2.12190386     -12.06174102       2.48123664
        H       2.68569500      -8.09238646       4.14547290
        H       2.43329579      -7.91168000      -1.56456659
        H       2.26368032     -11.97103749      -0.03692581
        H      -0.30289388     -11.03081473       2.99676260
        H       0.05404650      -8.59747583       4.07520644
        H      -0.14620660      -8.61151792      -1.35464113
        H      -0.25838864     -11.10457655      -0.36003217
        H      -0.85662057      -7.76723654       1.38227139
        H       0.26321351      -6.74399730       0.47215750
        H       0.31659006      -6.75672110       2.23713123
        C       0.17801499      -7.39347280       1.35339060
        H       4.83158681      -7.79302616       0.33628574
        H       2.98632816      -6.14645123       0.41958259
        C       4.37932750      -7.48955975       1.27484960
        C       3.37901781      -6.61143075       1.32052856
        H       3.01865674      -6.20459134       2.26225286
        H       4.86414560      -7.85327440       2.17514370
    End
End"""
)
mol.guess_bonds()
print(f"{mol.charge=}")
view(mol, direction="along_y", show_atom_labels=True, atom_label_type="Name")


# ## DFTB Settings

from scm.plams import Settings, AMSJob


def get_dftb_settings():
    # Engine settings: GFN-1xTB with implicit Tolene solvation
    settings = Settings()
    settings.input.DFTB.Model = "GFN1-xTB"
    settings.input.DFTB.Solvation.Solvent = "Toluene"
    return settings


print(AMSJob(settings=get_dftb_settings()).get_input())


# ## Geometry optimization of reactant state


def get_reactants_job(mol):
    settings = get_dftb_settings()
    settings.input.ams.Task = "SinglePoint"
    settings.input.ams.Properties.NormalModes = "Yes"
    job = AMSJob(molecule=mol, settings=settings, name="reactants_frequencies")
    return job


reactants_job = get_reactants_job(mol)
reactants_job.run()


# ## Set up PES Scan to find approximate transition state


def get_pesscan_job(mol):
    # Engine settings: GFN-1xTB with implicit Tolene solvation
    settings = get_dftb_settings()

    # AMS settings for PES scan (see online tutorial for choice of scan coordinates)
    settings.input.ams.Task = "PESScan"
    settings.input.ams.PESScan.ScanCoordinate.nPoints = "20"
    settings.input.ams.PESScan.ScanCoordinate.Distance = ["1 28 3.205 2.3", "29 25 3.295 1.5"]

    # Loosened convergence criteria
    settings.input.ams.GeometryOptimization.Convergence.Energy = "5.0e-5 [Hartree]"
    settings.input.ams.GeometryOptimization.Convergence.Gradients = "5.0e-3"
    settings.input.ams.GeometryOptimization.Convergence.Step = "5.0e-3 [Angstrom]"

    # setting up the AMS job with the molecule and settings object
    job = AMSJob(molecule=mol, settings=settings, name="pes_scan")
    return job

    # run the calulation


pesscan_job = get_pesscan_job(reactants_job.results.get_main_system())
pesscan_job.run()


# ## PESScan results

import matplotlib.pyplot as plt
from pprint import pprint

pesscan_results = pesscan_job.results.get_pesscan_results()
pprint(pesscan_results, depth=1)


fix, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].plot(list(range(pesscan_results["nPESPoints"])), pesscan_results["PES"])
ax[0].set_ylabel("Energy (hartree)")
ax[0].set_xlabel("Frame")
for i in range(pesscan_results["nRaveledScanCoords"]):
    ax[1].plot(
        list(range(pesscan_results["nPESPoints"])),
        pesscan_results["RaveledPESCoords"][i],
        label=pesscan_results["RaveledScanCoords"][i],
    )
ax[1].set_ylabel("bohr")
ax[1].set_xlabel("Frame")
ax[1].legend()


import numpy as np

highest_energy = np.max(pesscan_results["PES"])
highest_energy_index = np.argmax(pesscan_results["PES"])

# res['Molecules'] contains all the converged geometries (as PLAMS Molecules)
highest_energy_geo = pesscan_results["Molecules"][highest_energy_index]

print(f"Highest energy: {highest_energy:.6f} hartree, for geometry #{highest_energy_index+1}")
print("We will use that structure as starting point for TS search")
highest_energy_geo.guess_bonds()
view(highest_energy_geo, direction="along_y")


# ## Transition state search


def get_ts_job(mol):
    settings = get_dftb_settings()

    # AMS settings for TS search
    settings.input.ams.Task = "TransitionStateSearch"
    settings.input.ams.Properties.NormalModes = "Yes"
    settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"

    # setting up the AMS job with the molecule and settings object
    job = AMSJob(molecule=mol, settings=settings, name="ts_search")
    return job


ts_job = get_ts_job(highest_energy_geo)
ts_job.run()
ts_mol = ts_job.results.get_main_system()
ts_mol.guess_bonds()
view(ts_mol, direction="along_y")


def get_gibbs_energy(job: AMSJob, unit="hartree"):
    from scm.libbase import Units

    gibbs_energy = job.results.readrkf("Thermodynamics", "Gibbs free Energy", file="engine")
    gibbs_energy *= Units.convert("hartree", unit, 1.0)
    return gibbs_energy


def print_ts_results(job):
    energy = job.results.get_energy(unit="kcal/mol")
    gibbs_energy = get_gibbs_energy(job, unit="kcal/mol")
    frequencies = job.results.get_frequencies(unit="cm^-1")
    character = job.results.readrkf("AMSResults", "PESPointCharacter", file="engine")

    # print results
    print("== Results ==")
    print("Character   : {}".format(character))
    print("Energy      : {:.3f} kcal/mol".format(energy))
    print("Gibbs Energy: {:.3f} kcal/mol".format(gibbs_energy))
    for freq in frequencies:
        if freq < 0:
            print("Imaginary frequency : {:.3f} cm^-1".format(freq))
        else:
            break


print_ts_results(ts_job)


# ## Free energy barrier


def print_final_results(gibbs_energy_ts: float, gibbs_energy_reactants: float, T=298.15):
    barrier = gibbs_energy_ts - gibbs_energy_reactants  # kJ/mol
    NA = 6.022e23  # Avogadro's constant
    k_B = 1.380649e-23  # Boltzmann constant, J / (mol*K)
    h = 6.626e-34  # Planck's constant, J * s
    rate_constant = (k_B * T / h) * np.exp(-barrier * 1000 / (NA * k_B * T))
    print("Results summary")
    print("Reactants Gibbs Energy: {:.3f} kJ/mol".format(gibbs_energy_reactants))
    print("TS Gibbs Energy:        {:.3f} kJ/mol".format(gibbs_energy_ts))
    print("Free energy barrier:    {:.3f} kJ/mol".format(barrier))
    print("Rate constant:          {:.3f} s⁻¹ at {} K".format(rate_constant, T))


gibbs_energy_ts = get_gibbs_energy(ts_job, "kJ/mol")
gibbs_energy_reactants = get_gibbs_energy(reactants_job, "kJ/mol")
print_final_results(gibbs_energy_ts, gibbs_energy_reactants, T=298.15)