index.qmd

---
title: |
  **KIMMDY**\
  Reactive Molecular Dynamics\
  using\
  Kinetic Monte Carlo
# subtitle: Designing Software for Humans
date: today
date-format: June 3, 2024
author: Jannik Buhr, MBM
institute: Heidelberg Institute for Theoretical Studies
format:
  hits-revealjs:
    width: 1920
    height: 1080
    margin: 0.04
    code-line-numbers: false
    title-slide-attributes:
      data-background-image: ./www/helix-wide.png
      data-background-color: white
      data-background-size: contain
      data-background-position: center bottom -10%
    citations-hover: true
    reference-links: false
    preview-links: false
    progress: true
    slide-number: true
    chalkboard:
      buttons: true
      src: ./www/chalkboard.json
    multiplex: false
    css: ./www/styles.css
    footer: |
      <https://github.com/hits-mbm-dev/2024-hits-kimmdy>
footnotes-hover: true
reference-location: document
favicon: ./www/hits-logo.png
bibliography: ./references.bib
mermaid-format: png
nocite:
  - "@allaireQuarto2022c"
engine: jupyter
execute:
  eval: true
  echo: false
  warning: false
  cache: false
---

```{python}
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt

plt.rcParams['figure.dpi'] = 100
```

## Collagen has interesting chemistry under mechanical stress {background-image='./img/collagen.png' background-color=black}

## How do we simulating (large) molecular systems?

::: {.columns}
::: {.column width="50%"}
![](img/monkey.png){width="90%"}
:::
:::{.column width="50%"}
### Bonds do not exist.
:::{.r-stack .center}
![](./img/molecule-vdw.png)

![](./img/molecule-density.png){.fragment}

![](./img/molecule-cpk.png){.fragment}
:::

### But let's pretend they do. {.fragment}
:::
:::

## Molecular Dynamics / Molecular Mechanics -- Newton to the rescue!

:::: {.columns}
::: {.column width="40%"}
$$E_{total} = E_{bonded} + E_{nonbonded}$$
$$E_{bonded}=E_{\text{bond}}+E_{\text{angle}}+E_{\text{dihedral}}$$
$$E_{nonbonded} = E_{\text{electrostatic}} + E_{\text{van der Waals}}$$
:::
::: {.column width="30%"}
![](./img/molecule-cpk.png){#mol1}
:::
::: {.column width="30%"}
![](./img/molecule-cpk.png){#mol2}
:::
::::

### Popular Potential Energy Functions {.fragment}

<hr>

:::: {.columns}
::: {.column width="30%" .fragment}

:::{#oh}
![](./img/hyd-oh.png)
:::

```{python}
x = np.linspace(0.1, 10, 100)
d = 1
b = 1
y = d * (1 - np.exp(-b * (x - 1)))**2

plt.subplots(figsize=(3,2))
plt.plot(x, y, color='black', linewidth=5)
plt.box(False)
plt.xticks([0, b], ['0', '$r_0$'])
plt.ylabel('$E$')
plt.yticks([])
plt.show()
```

$$E_{morse}(r) = D \left(1 - e^{-\beta(r - r_0)}\right)^2$$
:::

::: {.column width="30%" .fragment}

:::{#water}
![](./img/water.png)
:::

```{python}
x = np.linspace(-2, 2, 100)
k = 0.1
y = 0.5 * k * x**2

plt.subplots(figsize=(2,2))
plt.box(False)
plt.xticks([0], ['$r_0$'])
plt.yticks([])
plt.plot(x, y, color='black', linewidth=5)
plt.show()
```

$$E_{harmonic}(r) = \frac{1}{2} k (r - r_0)^2$$

:::

::: {.column width="30%" .fragment}
```{python}
x = np.linspace(0, 2*pi, 100)
y = 1 + np.cos(2 * x - pi)

plt.subplots(figsize=(3,2))
plt.box(False)
plt.xticks([pi], [r'$\phi_0$'])
plt.yticks([])
plt.plot(x, y, color='black', linewidth=5)
plt.show()
```

$$E_{cosine}(\phi) = k_\phi (1 + \cos(n \phi - \phi_0))$$

:::
::::

:::{#meme-harmonic .fragment}
![](img/meme-harmonic.png){width="85%"}
:::

## Force Fields list the parameters for these potentials

e.g. AMBER99SB*-ILDNP [@cornellSecondGenerationForce1995; @alievMotionalTimescalePredictions2014]

::: {.columns}
::: {.column width="40%"}
```{.itp filename="forcefield.itp"}
#include "ffnonbonded.itp"
#include "ffbonded.itp"
```

```{.itp filename="ffnonbonded.itp"}
[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
Br          35      79.90    0.0000  A   0.00000e+00  0.00000e+00
C            6      12.01    0.0000  A   3.39967e-01  3.59824e-01
```

```{.itp filename="ffbonded.itp"}
[ bondtypes ]
; i    j  func       b0          kb
  C  C          1     0.1525   259408.0
  C  OS         1     0.1323   376560.0
[ angletypes ]
...
```

:::{#harmonic}
$$V_{harmonic}(r) = \frac{1}{2} k (r - r_0)^2$$
:::
:::
::: {.column width="50%"}
![](./img/molecule-cpk.png)
:::
:::

### But what about Chemistry? {.fragment .cent}

## Welcome to KIMMDY!

_Project with Kai Riedmiller and Eric Hartmann_

<https://hits-mbm-dev.github.io/kimmdy/>

{{< qrcode https://hits-mbm-dev.github.io/kimmdy/ qr-kmdy width=250 height=250 colorDark='#1E4287'  >}}

![](./img/ala-rad.png){#ala-rad}

<style>
#ala-rad {
  position: absolute;
  top: 50%;
  left: 0;
  width: 40%;
}
</style>

::: {.columns}

::: {.column width="20%"}
```bash
pip install kimmdy
```


```bash
edit kimmdy.yml
```
:::

::: {.column width="50%"}

```{dot}
//| fig-height: 8
//| fig-width: 3
digraph L {

  edge [minlen="1"]
  node [shape=rect fontname="Miriam Libre" style=filled fillcolor=white fontsize=8];

  md [ label="Run MD" fillcolor="#1e42870d"];
  reactions [ label="Query Reactions" fillcolor="#1e42873d"];
  mc [label="Choose Reaction" fillcolor="#1e42870d"];
  execute [label="Execute Reaction" fillcolor="#1e42870d"];

  md -> reactions -> mc -> execute -> md;

}
```
:::
::: {.column width="30%" .fragment}

```{python}
fig, ax = plt.subplots(figsize=(1,5))
# plt.axis('off')

HITS_BLUE = "#1E4287"
HITS_GREEN = "#019050"
HITS_MAGENTA = "#c3006b"
HITS_YELLOW = "#ffcc00"

counts = {
    '1': (0.5, HITS_BLUE),
    '3': (0.2, HITS_MAGENTA),
    '2': (0.1, HITS_YELLOW),
    '4': (0.2, HITS_GREEN),
}

bottom = np.zeros(1)

for k, (v,c) in counts.items():
    p = ax.bar(" ", v, width=1, label=k, bottom=bottom, color=c)
    bottom += v

    ax.bar_label(p, label_type='center', color='white')

plt.show()
```
:::
:::


:::{#dice .fragment}
![](img/dice.png){width="30%"}
:::

## Define your own workflow

::::{.columns}
:::{.column width="30%" .large}
```yaml
sequence:
  - equilibrium
  - mult: 4
    tasks:
    - reactions
    - equilibrium
```
:::

:::{.column width="70%"}
![](img/yml-warning.png){.fragment #yml-warning}
:::
::::

::::{.columns}
:::{.column width="30%" .large}
:::
:::{.column width="70%"}
:::
::::

## KIMMDY is user-friendly

::::{.columns}
:::{.column width="30%"}
- **checks** and **documentation**
  right in your editor
:::
:::{.column width="70%"}
:::{.r-stack}
![](img/yml-hover.png)

![](img/yml-completion.png){.fragment}
:::
:::
::::


## KIMMDY is extensible

::::{.columns}
:::{.column width="30%"}
### Builtin Reactions

- **Homolysis**
- **Hydrogen Atom Transfer** (HAT)

### Define your own reactions

- Simple python plugin
- Take **state** of the system
- Return **rates** and **recipes**
- KIMMDY takes care of topology details
:::
:::{.column width="70%" .fragment}
### Example: Peptide bond hydrolysis
![](./www/hydrolysis.mp4)
:::
::::


# Thank You!

:::{#otter}
![](./img/otter.png){width="30%"}
:::

## Thank You!

:::{#mbm-group}
![](./www/mbm-group.jpeg){width="80%"}
:::

:::{#qr}
Slides:
{{< qrcode https://github.com/hits-mbm-dev/2024-hits-kimmdy qr-slides width=256 height=256 colorDark='#1E4287' >}}
:::

::::{#bar}
:::profile
- [mbm.h-its.org](https://mbm.h-its.org){.bi-house}
- [hits-mbm-dev](hhttps://github.com/hits-mbm-dev){.bi-github}
- [jmbuhr](https://github.com/jmbuhr){.bi-github}
- [jmbuhr.de](https://jmbuhr.de){.bi-house}
:::
:::{#acknowledge}
This work was supported by the Klaus Tschira Foundation and has received funding from the European Research Council (ERC)
:::
::::

## References

:::{#refs}
:::

# Backup Slides

## Molecules

{{< mol-url ./examples/alanine_hat/just-ala.gro >}}

## Topology

::::{.column width=50%}
```{dot}
//| fig-height: 3
//| fig-width: 3
graph G {
  layout=neato
  overlap=true
  node [shape="circle"]
  "1 O"  -- "2 O" ;
}
```

```{dot}
//| fig-height: 3
//| fig-width: 3
graph G {
  layout=neato
  overlap=true
  node [shape="circle"]
  "1 H"  -- "2 O" ;
  "3 H"  -- "2 O" ;
}
```

```{dot}
//| fig-height: 3
//| fig-width: 5
graph G {
  layout=neato
  overlap=true
  rankdir=LR
  node [shape="circle"]
  "1"  -- "2" ;
  "2"  -- "3" [penwidth=5];
  "3"  -- "4" ;
}
```
::::