Molecular dynamics implementation cheat sheet

A cheat-sheet for implementing a molecular dynamics simulation
Source: "Nanosystems: Molecular Machinery, Manufacturing, and Computation"

The MM2 model (3.3.2.)

[math]\mathcal{V}_{s} = \frac{1}{2}k_{s}(r - r_{0})^{2}\lbrack 1 - k_{cubic}(r - r_{0})\rbrack[/math] … stretching (3.4)
[math]\mathcal{V}_{\theta} = \frac{1}{2}k_{\theta}(\theta - \theta_{0})^{2}\lbrack 1 + k_{sextic}(\theta - \theta_{0})^{4}\rbrack[/math] … bending (3.5)
[math]k_{s\bot} = k_{\theta}/r_{0}^{2}[/math] … per length angular stiffness [math]k_{s\bot} \approx k_{s}/20[/math] at [math](r_{0},\theta_{0})[/math] (3.6)
[math]\mathcal{V}_{\omega} = \frac{1}{2}\lbrack V_{1}(1 - \cos(1\omega)) + V_{2}(1 - \cos(2\omega)) + V_{3}(1 + \cos(3\omega))\rbrack[/math] … torsion (3.7)
[math]\mathcal{V}_{vdw} = \epsilon_{vdw}\lbrack 2.48 \times 10^{5}\exp( - 12.5\frac{r}{r_{vdw0}}) - 1.924(\frac{r}{r_{vdw0}})^{- 6}\rbrack[/math] … vdW attraction & non-bonded Pauli repulsion (3.8)
[math]\mathcal{V}_{s\theta} = k_{s\theta}(\theta - \theta_{0})\lbrack(r_{A} - r_{A0}) + (r_{B} - r_{B0})\rbrack[/math] … stretch-bend interaction (3.9)

---

[math]k_{s,C - C} = 440N/m[/math] — [math]r_{0,C - C} = 111.3pm[/math] … (Table 3.2.)
[math]k_{\theta,C - C - C} = 450zJ/rad^{2}[/math] — [math]\theta_{0,C - C - C} = 1.911rad = 109.47{^\circ}[/math] … (Table 3.3.)
[math]\epsilon_{vdw,Csp^{3}} = 357yJ[/math] — [math]r_{vdw,Csp^{3}} = 190pm[/math] … well depth [math]\epsilon_{vdw}[/math] varies widely for other atoms (~10x) (Table 3.1.)
[math]V_{1,C - C - C - C} = 1.39zJ,\ V_{2,C - C - C - C} = 1.88zJ,\ V_{3,C - C - C - C} = 0.65zJ[/math] … (Table 3.5.)
[math]k_{s\theta,C - C - C} = 1.2nN/rad[/math] … (Table 3.6.)
[math]k_{cubic} = ?m^{- 1}[/math] — [math]k_{sextic} = ?rad^{- 4}[/math]

[math]n … nano = 10^{-9}[/math] — [math]p … pico = 10^{-12}[/math]
[math]z … zepto = 10^{-21}[/math] — [math]y … yocto = 10^{-24}[/math]

Bonds under large loads (beyond MM2) (3.3.3.)

[math]\mathcal{V}_{morse} = D_{e}({1 - \exp\lbrack - \beta(r - r_{0})\rbrack}^{2} - 1)[/math] … tensile (3.10)
[math]\mathcal{V}_{lippincott} = D_{e}\lbrack 1 - \exp( - \frac{k_{s}r_{0}(r - r_{0})^{2}}{2D_{e}r})\rbrack,\ \ r \geq r_{0}[/math] … tensile – more accurate for larger distances (3.15)
For compressive loads [math]\mathcal{V}_{vdw}[/math] see (3.8) above [math]F_{vdw} = - \frac{\partial}{\partial r}\mathcal{V}_{vdw}[/math] … (3.16)
Note: Popular Lennard-Jones 6-12 potential is too steep in repulsive regime!
MM2 potential within 10% of experiment for [math]0.5r_{vdw}[/math] ([math]\gt 100yJ[/math] repulsion) 🙂

[math]\beta = \sqrt{k_{s}/(2D_{e})}[/math] … (3.13)
[math]D_{e} \approx D_{0} + \frac{\hslash}{2}\sqrt{k_{s}/\mu} = D_{0} + \frac{\hslash\omega}{2};\ \ \mu = \frac{m_{1}m_{2}}{m_{1} + m_{2}}[/math] … (3.14)
[math]D_{0}[/math] … potential well depth
[math]D_{e,C - C} = 556yJ[/math] — [math]k_{s,C - C} = 440N/m[/math] — [math]r_{0,C - C} = 152.3pm[/math] … (Table 3.8)

Verlet Integration (not in Nanosystems)

[math]\overset{\rightarrow}{x}(t + \Delta t) = 2\overset{\rightarrow}{x}(t) - \overset{\rightarrow}{x}(t - \Delta t) + \overset{\rightarrow}{a}(t)\Delta t^{2} + O(\Delta t^{4})[/math]
[math]\overset{\rightarrow}{a}(t) = m\nabla_{\overset{\rightarrow}{x}}\mathcal{V}(t)[/math]

MM2 vs MM3 (Chapter 3.3.2.g.)

★ MM2 prioritizes accuracy in energy and geometry
★ MM3 prioritizes accuracy of vibration frequencies
★ MM3 predicts greater angle-bending stiffness by ~1.5x or more
★ MM3 has a more complex functional form – e.g. it adds:
– stretch-torsion interaction, cubic bending, quartic stretching
★ MM3 has lower energies forces and stiffness in deep repulsive regime by ~10%
★ MM2 overall seems like a more conservative (safe side wrong) choice for diamondoid nanomachinery

Bond cleavage & radical coupling (3.4.2)

The Morse potential [math]\mathcal{V}_{morse}[/math] can approximate homolytic reactions
[math]\mathcal{V}_{anti - morse} = \frac{1}{2}D_{e}({1 + \exp\lbrack - \beta(r - r_{0})\rbrack}^{2} - 1)[/math] … unpaired spins (3.21)
But unpaired spins being the limiting factor in reaction rate should
typically be avoidable in piezomechanosynthesis.
See Nanosystems Chapter 8.4.3.b. Radical coupling and intersystem crossing.

Nonbonded & large compression limits (3.3.3.b.)

[math]k_{s,vdw} \approx \frac{12.5}{r_{vdw0}}F_{vdw} \approx 3.5 \times 10^{10}m^{- 1} \cdot F_{vdw}[/math] … (3.18)
[math]\mathcal{V}_{vdw} \approx 0.08r_{vdw0}F_{vdw} \approx 2.9 \times 10^{- 11}m \cdot F_{vdw}[/math] …(3.19)

Abstraction reactions (Chapter 3.4.3.)

[math]\mathcal{V}_{LEPS}[/math] extended London-Eyring-Polyano-Sato potential – details omitted here

Continuum models of van der Waals attraction (Chapter 3.5.1.)

Hamaker constant – details omitted here